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The  finite  element  method  may  relieve  some  of  the 
problems  associated  with  numerical  modelling  of  the  atmos¬ 
phere.  The  method  reduces  the  problem  of  nonlinear  computa¬ 
tional  instability,  allows  arbitrary  placement  of  grid  points, 
and  offers  greater  flexibility  in  the  handling  of  boundary 
conditions.  In  addition,  the  finite  element  method  has  been 
shown  to  be  more  accurate  for  some  types  of  problems. 

This  research  concentrates  on  the  ability  of  the  fin¬ 
ite  element  method  to  serve  as  the  means  of  solving  the  equa¬ 
tions  which  describe  flow  in  the  atmosphere  and  provides 
answers  as  to  the  type  of  finite  element  approximation  best 
suited  for  meteorological  research.  Five  different  finite 
elements  using  both  linear  and  quadratic  interpolation  are 
tested  in  the  space  domain.  Several  time  differencing 
schemes  are  also  tested  with  the  elements  in  order  to  deter¬ 
mine  the  most  accurate  and  efficient  configuration.  In  add¬ 
ition,  the  concepts  of  lumped  and  consistent  mass  are  examined 
and  tested  as  well  as  alternative  methods  of  handling  the 
computer  implementation  of  the  method. 

It  is  found  that  the  finite  element  method  using  the 
four-node  bilinear  rectangular  element  provides  the  most  ac- 
cruate  handling  of  the  spatial  derivatives.  Combined  with 
the  Crank-Nicholson  time  difference  scheme  using  consistent 
mass,  the  four-node  element  is  not  only  the  most  accurate  in 
the  handling  of  advective  flow  but  also  controls  the  growth 
of  gravity  waves  the  longest.  This  study  also  shows  that 
when  explicit  time  difference  schemes  such  as  the  Leap-frog 
scheme  are  used  with  consistent  mass  the  degradation  in  re¬ 
sults  when  compared  to  the  Crank-Nicholson  method  is  minimal. 

From  the  standpoint  of  computational  efficiency,  the 
finite  element  method  is  shown  to  be  competitive  with  certain 
finite  difference  methods  under  certain  special  conditions. 

For  constant  grid  spacing  the  finite  element  version  of  a 
simple  vorticity-stream  function  model  using  the  four-node 
bilinear  element  can  be  converted  into  an  expression  identical 
to  the  Arakawa  Jacobian  in  space.  However,  unlike  Arakawa's 
method,  the  finite  element  method  includes  consistent  mass, 


a  more  accurate  distribution  of  the  mass  than  found  in  con¬ 
ventional  finite  difference  methods  which  use  the  lumped 
mass  concept.  This  special  finite  difference  form  of  the 
finite  element  method  requires  slightly  more  computation  but 
is  more  accurate  than  its  finite  difference  counterpart. 

It  is  evident  that  further  study  of  the  finite  ele¬ 
ment  method  as  an  alternative  to  finite  difference  methods 
is  warranted.  Its  use  in  the  future  may  solve  many  existing 
problems  in  current  atmospheric  models  as  well  as  providing 
more  accuracy  in  the  solutions  of  the  equations. 
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ABSTRACT 


\  • 

U 

The  finite  element  method  may  relieve  some  of  the 
problems  associated  with  numerical  modelling  of  the  atmos¬ 
phere.  The  method  reduces  the  problem  of  nonlinear  computa¬ 
tional  instability,  allows  arbitrary  placement  of  grid  points, 
and  offers  greater  flexibility  in  the  handling  of  boundary 
conditions.  In  addition,  the  finite  element  method  has  been 
shown  to  be  more  accurate  for  some  types  of  problems. 

This  research  concentrates  on  the  ability  of  the 
finite  element  method  to  serve  as  the  means  of  solving  the 
equations  which  describe  flow  in  the  atmosphere  and  provides 
answers  as  to  the  type  of  finite  element  approximation  best 
suited  for  meteorological  research.  Five  different  finite 
elements  using  both  linear  and  quadratic  interpolation  are 
tested  in  the  space  domain.  Several  time  differencing  schemes 
are  also  tested  with  the  elements  in  order  to  determine  the 
most  accurate  and  efficient  configuration.  In  addition,  the 
concepts  of  lumped  and  consistent  mass  are  examined  and  tested 
as  well  as  alternative  methods  of  handling  the  computer  im¬ 
plementation  of  the  method. 

It  is  found  that  the  finite  element  method  using  the 
four-node  bilinear  rectangular  element  provides  the  most 


accurate  handling  of  the  spatial  derivatives.  Combined  with 
the  Crank-Nicholson  time  difference  scheme  using  consistent 
mass,  the  four-node  element  is  not  only  the  most  accurate  in 
the  handling  of  advective  flow  but  also  controls  the  growth 
of  gravity  waves  the  longest.  This  study  also  shows  that 
when  explicit  time  difference  schemes  such  as  the  Leap-frog 
scheme  are  used  with  consistent  mass  the  degradation  in 
results  when  compared  to  the  Crank-Nicholson  method  is  min¬ 
imal  . 

From  the  standpoint  of  computational  efficiency,  the 
finite  element  method  is  shown  to  be  competitive  with  certain 
finite  difference  methods  under  certain  special  conditions. 

For  constant  grid  spacing  the  finite  element  version  of  a 
simple  vorticity- stream  function  model  using  the  four-node 
bilinear  element  can  be  converted  into  an  expression  identical 
to  the  Arakawa  Jacobian  in  space.  However,  unlike  Arakawa's 
method,  the  finite  element  method  includes  consistent  mass, 
a  more  accurate  distribution  of  the  mass  than  found  in  conven¬ 
tional  finite  difference  methods  which  use  the  lumped  mass 
concept.  This  special  finite  difference  form  of  the  finite 
element  method  requires  slightly  more  computation  but  is 
more  accurate  than  its  finite  difference  counterpart. 

It  is  evident  that  further  study  of  the  finite  ele¬ 
ment  method  as  an  alternative  to  finite  difference  methods 
is  warranted.  Its  use  in  the  future  may  solve  many  existing 
problems  in  current  atmospheric  models  as  well  as  providing 
more  accuracy  in  the  solutions  of  the  equations. 
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THE  USE  OF  THE  FINITE  ELEMENT  METHOD 
IN  METEOROLOGICAL  MODELLING 

CHAPTER  I 

INTRODUCTION 

A  large  percentage  of  current  meteorological  research 
involves  the  development  of  numerical  models.  The  majority 
of  those  models  are  designed  using  the  finite  difference 
method  (FDM)  to  represent  the  derivatives  in  the  differential 
equations.  The  use  of  the  FDM  in  meteorological  models  leads 
to  difficulties  in  several  areas.  For  example,  many  differ¬ 
encing  schemes  result  in  damping  of  important  waves,  inaccur¬ 
acies  in  phase  speeds,  or  computational  instability  (Messinger 
&  Arakawa,  1976).  In  addition,  finite  difference  methods  are 
not  equipped  to  handle  irregular  boundaries.  Finally,  one  of 
the  most  difficult  situations  for  the  FDM  is  grid  nesting,  an 
important  feature  in  numerical  forecasting  which  allows  for 
grid  refinement  in  areas  of  interest  so  that  smaller  scale 
features  can  be  identified  and  predicted  with  greater  accuracy. 

One  method  which  may  relieve  some  of  the  problems  of 
finite  differencing  is  the  finite  element  method  (FEM) .  The 
FEM  employs  local  interpolating  functions  over  finite  sub- 


domains  which  are  connected  together  to  approximate  the  glo¬ 
bal  system.  The  governing  partial  differential  equations 
are  expressed  in  integral  form.  Normally  this  is  accomplish¬ 
ed  using  the  Galerkin  method  which  allows  for  easy  derivation 
of  the  integral  form  of  complex  systems  of  partial  differen¬ 
tial  equations.  Details  on  the  Galerkin  method  are  given  in 
Appendix  A. 

Because  integral  equations  are  solved,  the  FEM  satis¬ 
fies  important  global  conservation  laws  associated  with  the 
model  equations,  regardless  of  the  shape  of  the  domain.  Thus, 
computational  instability  is  much  less  of  a  problem.  The 
nodal  (grid)  points  may  be  placed  arbitrarily  so  that  fine 
resolution  can  be  used  in  an  area  of  interest,  with  coarse 
resolution  in  areas  of  weak  gradients.  Higher  order  approx¬ 
imations  are  easily  used  without  special  boundary  conditions. 
There  is  greater  flexibility  in  the  handling  of  boundary  con¬ 
ditions.  Finally,  the  FEM  has  been  shown  to  be  more  accurate 
than  the  finite  difference  method  for  some  types  of  problems. 

In  summary,  the  finite  element  method  and  finite 
difference  method  are  fundamentally  different.  Finite  dif¬ 
ference  equations  are  derived  from  truncated  Taylor  Series 
expansions.  The  assumption  is  that  these  approximations 
offer  sufficient  accuracy  for  the  representation  of  differ¬ 
ential  equations  which  govern  physical  processes.  In  some 
finite  difference  methods,  the  concentration  is  on  the  partial 
differential  equation  and  the  approximation  of  derivatives  at 


a  point.  The  finite  element  method  is  global  in  nature;  it 
is  a  variational  method  which  uses  integral  equations. 

This  research  concentrates  on  the  ability  of  the  FEM 
to  serve  as  the  means  for  solving  the  equations  which  des¬ 
cribe  flow  in  the  atmosphere  and  provides  answers  as  to  the 
type  of  finite  element  approximation  best  suited  for  meteor¬ 
ological  research.  The  FEM  is  examined  in  several  different 
ways  in  order  to  determine  the  more  accurate  and  cost-effec¬ 
tive  ways  in  which  it  can  be  employed.  Five  different  ele¬ 
ments  using  both  linear  and  quadratic  interpolation  are  test¬ 
ed  in  the  space  domain.  These  elements  and  interpolation 
functions  are  discussed  in  Appendix  A.  Several  time  differ¬ 
encing  schemes  are  tested  with  the  elements  in  order  to  det¬ 
ermine  the  most  accurate  and  efficient  configuration.  In 
addition,  the  effect  of  using  consistent  versus  lumped  mass 
is  tested.  This  difference  is  concerned  with  the  handling 
of  the  time  derivative.  When  lumped  mass  is  used,  the  mass 
of  an  element  is  said  to  be  equally  distributed  to  the  nodes 
(or  mesh  points) .  In  contrast,  consistent  mass  means  that 
there  is  a  weighted  distribution  of  the  mass.  Lumped  mass 
is  analogous  to  the  finite  difference  method  and  is  a  com¬ 
promise  which  can  be  made  in  the  FEM.  Consistent  mass  natur¬ 
ally  arises  in  the  finite  element  formulation.  This  concept 
is  discussed  in  detail  in  Chapter  III. 

The  tests  described  above  are  performed  using  three 
different  formulations  for  flow  in  a  periodic  channel.  First 
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a  vorticity-stream  function  formulation  with  known  analytic 
solution  is  used  for  all  the  tests  described  above.  Then  the 
same  type  of  flow  is  studied  using  the  penalty  method  for 
comparison.  In  the  penalty  method  the  pressure  terms  are 
eliminated  in  the  functional;  however,  the  model  allows  for 
approximate  pressure  calculations.  Only  the  tests  involving 
element  accuracy  are  repeated  for  the  penalty  model.  Finally, 
the  shallow  water  equations  are  used  to  study  element  accur¬ 
acy,  time  differencing  methods,  and  the  conservative  proper¬ 
ties  of  the  FEM. 

Most  of  the  published  research  on  the  FEM  by  meteor¬ 
ologists  has  involved  a  specific  problem  with  the  FEM  applied 
in  a  specific  manner.  There  has  been  no  careful  investiga¬ 
tion  as  to  what  type  of  finite  element  approximation  is  best 
suited  to  meteorological  problems.  This  research  provides  a 
strong  foundation  for  further  application  of  the  FEM  in  ad¬ 
vanced  meteorological  models.  The  results  presented  here 
show  that  it  is  feasible  to  apply  the  FEM  to  non-linear  flow 
problems  where  advection  is  significant.  A  four-node  quadri¬ 
lateral  element  with  linear  interpolation  provides  the  best 
results  and  the  FEM  is  shown  to  be  more  accurate  than  the  FDM. 

Chapter  II  provides  insight  into  the  historical  devel¬ 
opments  of  the  FEM,  as  well  as  the  relationship  of  this  work 
to  published  research.  Chapter  III  gives  the  methodology 
including  a  discussion  of  the  governing  equations  and  numer¬ 
ical  procedures.  The  results  are  presented  in  Chapter  IV  and 


conclusions  and  remarks  in  Chapter  V.  There  are  appendices 
describing  the  finite  element  approximation  in  detail  and 
the  numerical  stability  for  the  four-node  rectangle. 


CHAPTER  II 


BACKGROUND 

Historical  Development 

The  variational  approach  to  the  solution  of  differ¬ 
ential  equations  is  the  basis  of  such  discoveries  as  Hamil¬ 
ton's  Principle  on  the  mechanics  of  moving  particles,  the 
Schrodinger  wave  equation,  and  contributed  to  Einstein's  dev- 
elopmentof  the  theory  of  relativity  (Simmons,  1972).  The 
variational  method  of  solving  differential  equations  is  quite 
elegant  for  certain  special  cases;  however,  its  general  appli¬ 
cability  has  been  delayed  until  modern  times  due  to  technical 
problems.  Approximate  methods  were  introduced  by  Rayleigh 
and  Ritz  during  the  late  1800 's  in  which  the  variables  were 
expressed  as  a  linear  combination  of  some  approximating  func¬ 
tions  and  the  integral  over  the  domain  of  the  error  was  mini¬ 
mized.  The  method  received  considerable  attention  but  it  is 
difficult  to  apply  for  problems  with  irregular  boundaries  or 
complex  boundary  conditions  and  cannot  be  used  for  nonlinear 
problems . 

A  break  occurred  in  the  1940 's  when  Courant  solved 
the  St.  Vincent  torsion  problem  by  approximating  the  warping 


function  through  a  combination  of  linear  triangles  assembled 
over  the  domain  and  formulating  the  problem  using  the  prin¬ 
ciple  of  minimum  potential  energy  (Oden,  1972).  It  is  found 
that  if  the  elements  are  relatively  small,  the  variation  of 
functions  within  them  can  be  adequately  represented  by  a 
low-order  polynomial.  If  certain  continuity  requirements 
were  met  along  the  element  boundaries ,  complex  systems  could 
be  reasonably  solved. 

The  method  came  to  the  forefront  in  the  1960 's  when 
its  use  became  widespread,  particularly  in  areas  such  as 
airframe  design.  Since  the  advent  of  large  scale  computer 
systems,  the  use  of  the  FEM  has  grown  rapidly  in  several 
areas  of  engineering;  however,  it  has  only  been  in  the  1970's 
that  interest  in  the  method  has  begun  to  develop  among  meteor- 
olgists. 


Current  Literature 


Other  than  the  University  of  Oklahoma,  School  of 
Meteorology,  there  are  four  main  groups  investigating  the 
application  of  the  FEM  to  meteorological  problems. 

Naval  Postgraduate  School,  Monterey,  California 
Meteorological  Office,  Bracknell,  England 
Atmospheric  Environment  Service,  Quebec 
Lawrence  Livermore  Laboratory,  Livermore,  California 
Additionally,  there  has  been  some  related  research  by  civil 
engineers  and  oceanographers  such  as  Kawahara  (1977  a,  b,  c) , 
Platzman  (1978) ,  and  Fix  (1975)  who  have  employed  the  shallow 
water  equations.  This  survey  of  current  literature  will 
provide  a  capsule  description  of  the  state  of  meteorological 
FEM  research  and  emphasize  the  need  for  this  study. 

Work  at  the  Naval  Postgraduate  School  culminated  with 
the  publication  of  Kelly  and  Williams  (1976)  which  was  a  study 
of  barotropic  flow  in  a  periodic  channel.  The  model  used  was 
the  shallow  water  equations.  The  two  space  dimensions  were 
approximated  using  the  linear  FEM  over  triangular  elements. 

The  time  integration  was  performed  using  centered  (Leap-frog) 
finite  differencing  with  consistent  mass.  Several  grid 
arrangements  were  attempted  using  well-behaved  initial  con¬ 
ditions  (wave  number  1) .  The  results  displayed  after  48  hours 
of  time  integration  were  disappointing.  There  was  mild  im¬ 
provement  only  with  the  addition  of  a  diffusion  term  or  when 
a  fine  mesh  was  employed.  They  tested  two  different  arrange- 


ments  of  the  diagonal  slant  of  the  triangular  elements  but 
did  not  identify  whether  the  noise  problems  were  associated 
with  space  or  time  differencing. 

The  FEM  research  at  the  Meteorological  Office,  Brack¬ 
nell,  England  has  been  published  by  Cullen.  His  first  work 
(Cullen,  1973)  showed  that  the  rectangular  finite  element 
with  bilinear  interpolation  and  consistent  mass  leap-frog 
time  differencing  results  in  a  more  restrictive  (by  a  factor 
of  /3)  maximum  time  step  than  the  one  required  when  centered 
finite  differencinq  is  used.  Additionally;  he  found  that  a 
16  x  16  finite  element  grid  gave  results  comparable  to  a  32 
x  32  finite  difference  grid  with  second  order  differencing. 
His  next  publication  (Cullen,  1974a)  described  his  investiga¬ 
tion  of  the  shallow  water  equations  using  the  same  model 
problem  as  used  by  Grammeltvedt  (1969)  .  Cullen  found  that 
linear  FEM  approximation  and  triangular  elements  gave  super¬ 
ior  results  to  the  finite  difference  models  tested  by  Gramm¬ 
eltvedt.  Specifically,  the  FEM  handled  wave  number  1  to  5 
(16  grid  points  along  the  channel)  with  90  percent  accuracy 
or  better  while  the  finite  difference  models  with  double  the 
resolution  only  achieved  90  percent  accuracy  for  wave  numbers 
1  to  3.  Additionally,  he  found  that  his  FEM  model  generated 
extraneous  waves  of  two  grid  interval  length. 

Cullen  (1974b)  extended  his  work  with  the  shallow 
water  equations.  The  equations  were  expressed  in  spherical 
coordinates  and  were  solved  over  a  linear  triangular  mesh 
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where  the  globe  was  divided  into  large  icosahedrons  which 
were  further  subdivided  into  many  small  triangles.  The  re¬ 
sulting  global  grid  had  1002  nodes  (8-10  degrees  apart) . 
Initial  conditions  were  the  same  as  used  by  Phillips  (1959) 
but  because  the  icosahedral  grid  is  unsymmetric,  Cullen  was 
forced  to  solve  the  system  over  the  entire  globe  rather  than 
over  one  octant  as  Phillips  (1959)  did.  To  discretize  the 
time  domain,  Cullen  used  the  Leap-frog  method.  Rather  than 
using  a  pure  FEM,  Cullen  added  a  17-point  spatial  filter  dev¬ 
eloped  by  Shapiro  (1971).  With  this  model,  Cullen  showed 
that  the  1002-node  finite  element  model  performed  better  than 
a  finite  difference  model  using  4032  points.  However,  he  did 
encounter  noise  problems  in  the  model,  especially  at  the  in¬ 
tersections  of  the  icosahedrons.  He  attributed  the  noise  to 
the  spatial  discretization. 

Because  of  the  noise  problems,  Cullen  (1976)  contin¬ 
ued  his  research  by  looking  at  the  shallow  water  model  in 
combination  with  artificial  smoothing  methods.  Four  smooth¬ 
ing  schemes  were  tried,  including  fourth  order  nonlinear 
diffusion  terms  in  all  three  equations,  the  Shapiro  filter 
used  in  the  previous  paper  (Cullen,  1976b) ,  the  addition  of 
nonlinear  diffusion  terms  in  the  u  and  v  equations  only  and 
the  Sadourney  (1973)  method.  The  latter  method  provided  the 
best  results.  It  is  designed  to  damp  the  gravity  waves. 

In  summary,  it  is  clear  that  Cullen  made  important 
advances  in  the  use  of  the  FEM  but  it  should  be  noted  that 
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virtually  all  of  his  work  was  done  using  the  linear  triangu¬ 
lar  element  with  Leap-frog  time  differencing. 

Some  of  the  more  advanced  uses  of  the  FEM  are  being 
made  by  the  Atmospheric  Environment  Service  of  Canada.  The 
principal  researchers  have  been  Staniforth,  Mitchell,  and 
Daley.  Apparently,  their  original  intent  was  to  develop  an 
operational  regional  baroclinic  model  using  the  FEM.  Stani¬ 
forth  and  Mitchell  (1977)  studied  the  use  of  semi-implicit 
time  integration  using  the  scheme  of  Kwizak  and  Robert  (1971) . 
That  scheme  was  used  in  conjunction  with  a  barotropic  model 
where  both  bilinear  and  biquadratic  interpolation  functions 
were  employed.  Compared  to  the  finite  difference  results  of 
Kwizak  and  Robert,  the  FEM  provided  superior  results  espec¬ 
ially  in  terms  of  noise  reduction  and  reduced  damping.  The 
better  FEM  result  was  given  by  the  quadratic  interpolation 
functions;  however,  their  further  use  was  ruled  out  by  the 
researchers  because  of  the  increased  computational  costs. 

Staniforth  and  Daley  (1977)  expanded  to  a  three-dim¬ 
ensional,  primitive  equation  model  where  only  the  vertical 
coordinate  was  discretized  using  the  FEM.  The  horizontal 
domain  was  approximated  through  an  existing  spectral  model. 
Staniforth  and  Mitchell  (1978)  returned  to  the  shallow  water 
equations  to  investigate  the  effect  of  variable-resolution 
grids.  They  demonstrated  that  the  forecast  error  is  signif¬ 
icantly  reduced  when  a  smoothly-varying  mesh  size  is  used  to 
refine  the  grid  rather  than  when  the  grid  size  changes 
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through  a  discrete  jump.  Within  an  area  of  interest  the  mesh 
size  was  held  at  a  constant  fine  resolution  and  became  more 
and  more  coarse  away  from  that  area.  Again,  for  this  problem 
they  used  the  semi-implicit  time  integration  scheme  and  bi¬ 
linear  interpolation  over  rectangles.  An  interesting  side¬ 
light  to  this  research  is  their  use  of  Simpson's  rule  for  the 
element  integrations  rather  than  the  more  widely  used  Gauss 
quadrature  method.  Simpson's  rule  slightly  compromises 
accuracy  for  faster  speed. 

The  most  recent  paper  by  the  Canadian  group  (Stani- 
forth  and  Daley,  1979)  describes  a  limited  area,  three-dimen¬ 
sional,  baroclinic  finite  element  model.  The  bilinear  inter¬ 
polation  previously  used  was  generalized  to  handle  the  three 
dimensional  model.  The  model  resolution  was  relatively 
coarse  with  seven  levels  in  the  vertical  and  285  km  between 
nodes  in  the  fine  mesh  area.  The  results  were  encouraging, 
showing  this  model  to  be  competitive  with  an  operational  29 
wave  spectral  model  for  forecast  periods  up  to  48  hours. 

In  summary,  this  Canadian  group  has  made  advanced 
application  of  the  FEM  using  the  Kwizak  and  Robert  semi-im- 
plicit  time  integration  scheme  with  the  bilinear  finite  ele¬ 
ment  interpolation  functions. 

In  contrast  to  the  Canadian  work  on  large  scale 
models,  the  goal  of  researchers  at  the  Lawrence  Livermore 
Laboratory  is  to  construct  a  three  dimensional  boundary  layer 
model  based  on  the  FEM.  Gresho,  Lee,  and  Sani  (1977)  studied 
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the  tradeoffs  between  lumped  and  consistent  mass  for  linear 
advection  and  diffusion  problems.  They  found  that  when  lump¬ 
ed  mass  is  employed,  the  phase  speed  and  amplitude  exhibit 
much  greater  errors.  They  concluded  that  lumped  mass  "ser¬ 
iously  compromises  the  accuracy  of  the  FEM"  and  recommended 
that  the  effect  of  lumped  mass  be  investigated  for  the  non¬ 
linear  equations. 

Huyakorn,  et^  al.  (1978)  compared  four  elements  (six- 
node  quadratic  triangle  as  well  as  the  four-,  eight-,  and 
nine-node  rectangles)  using  the  mixed-order  interpolation 
method  to  investigate  steady  flow  through  a  sudden  expansion 
and  steady  free  thermal  convection  in  a  square  cavity.  Their 
results  showed  that  for  those  types  of  flow,  the  nine-node 
Lagrangian  element  gave  the  best  accuracy.  They  found  that 
the  accuracy  when  using  the  six-node  quadratic  triangle  was 
highly  dependent  on  the  arrangement  of  the  triangles.  The 
least  accuracy  was  given  by  the  eight-node  serendipity  ele¬ 
ment.  (see  Appendix  A  for  a  discussion  of  the  elements) . 

This  result  was  used  by  Gresho,  et  al.  (1978)  to  test  a  pre¬ 
dictor-corrector  method  for  time  integration.  The  method 
consisted  of  the  Adams-Bashforth  scheme  as  the  predictor  and 
the  trapezoidal  rule  as  the  corrector.  The  unique  feature  of 
their  algorithm  is  the  provision  for  calculation  of  the  time 
truncation  error  and  the  subsequent  adjustment  of  the  time 
step  to  reduce  the  error.  This  experiment  was  conducted  using 
mixed-order  interpolation  for  flow  starting  from  rest  in  a 
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channel  with  a  sudden  expansion.  Lumped  mass  was  again  tried 
with  much  better  results;  however,  this  improvement  was 
thought  to  be  due  to  the  much  lower  Reynold's  number  of  the 
latter  experiment. 

Gresho,  et  al_.  (1978)  summarize  much  of  their  earlier 
work  on  the  Navier-Stokes  equations  for  flow  in  a  channel 
with  a  cylindrical  obstruction.  The  major  thrust  of  this 
work  was  to  demonstrate  the  usefulness  of  their  semi-implicit 
time  differencing  scheme  described  earlier. 

Other  related  research  accomplished  by  oceanographers 
includes  a  work  by  Fix  (1975)  dealing  with  mesoscale  ocean 
flows.  He  described  a  simple  vorticity-stream  function  model 
using  quadratic  and  cubic  triangular  elements  in  space.  No 
computational  results  are  presented.  Similar  work  is  present¬ 
ed  by  Platzman  (1978)  who  used  a  linearized  set  of  primitive 
equations  for  a  coarse  resolution  ocean  circulation  problem. 
His  grid  was  identical  to  that  of  Cullen  (1974b) .  Several 
papers  dealing  with  tidal  flow  using  the  shallow  water  equa¬ 
tions  have  been  presented  by  Kawahara  (1977  a,  b,  c) .  Inter¬ 
ested  in  civil  engineering  aspects,  he  has  modelled  actual 
harbors,  tributaries,  and  lakes  using  the  FEM.  The  ease  with 
which  the  FEM  handles  irregular  boundaries  such  as  Tokyo 
Harbor  is  amply  demonstrated. 

The  research  described  above  helps  to  define  some  of 
the  large  problem  areas  requiring  investigation.  Other  than 
the  research  at  Lawrence  Livermore  Laboratory,  most  of  the 
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published  works  have  described  the  use  of  a  particular  fin¬ 
ite  element  (generally  the  linear  triangle)  with  one  time 
integration  scheme. 

At  the  University  of  Oklahoma,  School  of  Meteorology, 
research  was  begun  by  looking  closely  at  the  FEM  itself  be¬ 
fore  applying  the  method  to  a  sophisticated  model.  Using  a 
finite  difference  equivalent  to  one  of  the  FEM  models  tested 
here  as  well  as  some  FDM  models,  Sasaki  and  Reddy  (1977) 
studied  several  well-known  time  differencing  schemes  as  well 
as  the  variational  adjustment  technique  of  Sasaki  (1976) . 

They  did  not  specifically  consider  the  concept  of  lumped  and 
consistent  mass  but  performed  their  testing  using  the  FDM  as 
it  is  normally  formulated.  They  found  that  the  FDM  form  of 
the  bilinear  FEM  combined  with  the  Crank-Nicholson  time  dif¬ 
ferencing  scheme  provided  the  best  results.  This  scheme  is, 
in  fact,  a  consistent  mass  scheme;  the  others  tested  were  not. 
Sasaki  and  Chang  (1979)  carried  this  work  further  by  using  a 
consistent  mass  operator  with  some  of  the  finite  difference 
schemes  previously  tested.  They  found  improvement  in  the 
solution  for  each  of  the  schemes  when  consistent  mass  was 
employed . 

This  research  is  a  continuation  of  the  work  described 
above  but  is  broader  in  scope.  Here,  the  emphasis  is  on  the 
FEM  as  it  applies  to  the  equations  governing  large  scale  flow 
in  the  atmosphere.  Using  advection-dominated  flow  with  the 
nonlinear  equations,  the  discussion  which  follows  will  provide 
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insight  into  the  relative  accuracies  of  several  commonly  used 
finite  elements,  the  tradeoffs  associated  with  several  well- 
known  time  differencing  schemes,  the  relative  merits  of  lump¬ 
ed  and  consistent  mass,  and  computational  efficiency. 
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:  j  Governing  Equations  and  Finite  Element  Models 

[  The  objective  of  this  study  is  to  examine  the  FEM 

approximation  in  space  coupled  with  various  time  discretiza¬ 
tion  methods  in  order  to  determine  the  most  accurate  and 

{ 

\ 

t  efficient  application  of  the  FEM  for  use  in  modelling  the 

Accordingly,  three  dif- 
a  channel  will  be  used. 
The  first  two  models  are  derived  from  the  following  equation 
set  which  describes  inviscid,  incompressible  flow  in  a  chan¬ 
nel  with  no  Coriolis  force: 


large-scale  flow  of  the  atmospl{e^e. 
ferent  models  which  describe  flow  \n 
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Here  u  and  v  are  the  velocity  components  in  the  x  and  y  dir¬ 
ections  respectively,  p8  is  the  constant  density,  p  is  the 
pressure,  t  is  the  time,  and  (x,y)  denotes  a  point  in  ft,  an 
open  bounded  region  in  two-dimensional  Euclidean  space  with 
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boundary  denoted  by  3C2. 


Vorticity  -  Stream  Function  FEM  Model 


The  first  model  investigated  is  the  vorticity-stream 
function  model.  The  stream  function  (ij/)  is  defined  by 

dip  3^ 

u  ~  "  3y‘  v  =  3x  *  (3) 


The  vorticity  (C)  is  given  by 


3v  _  3u 
Tj?  cfy 


The  use  of  the  stream  function  automatically  satisfies  (2)  so 
that  (1),  (2),  (3),  and  (4)  may  be  combined  giving 

.jS.  +  =  0/ 


v2ip=  c. 


where  J(i i|>,?)  =  ^  is  the  Jacobian  operator. 

If  the  domain  8  is  subdivided  into  small  elements,  Q 

e 

the  vorticity  and  stream  function  may  be  approximated  over  an 
element  (e)  by 

ce  a  E?®  Ni  and  ipe  =  EipJ  Ni  .  (6) 

Here,  the  subscript  i  denotes  the  node  or  grid  point  i  and 
represents  the  element  interpolation  or  shape  functions. 
The  interpolation  functions  have  the  property  that  all  but 
one  are  zero  at  a  given  node;  the  one  corresponding  to  that 


node  is  unity.  Using  (6) ,  the  application  of  the  Galerkin 
method  to  (5)  results  in 
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where  M? .  =  ff  N.N .  dx  dy 
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It  should  be  noted  that  the  matrices  [M  ]  called  consistent 

0 

mass  matrix  and  [K  ]  are  constant  and  need  to  be  computed  and 

stored  only  once  for  each  element;  however,  the  nonlinear 
0 

matrix  [J  ]  must  be  recomputed  on  each  iteration  and/or  time 
step,  depending  on  the  time  integration  scheme  used.  Compu¬ 
tational  details  are  given  in  later  sections  of  this  chapter. 


Penalty  FEM 

An  alternative  to  the  vorticity  -  stream  function 
formulation  of  (1)  and  (2)  is  the  penalty  method,  which  is 
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also  known  as  the  weak-constraint  method  (Courant,  1943  and 
Reddy,  1979a,  b) .  In  the  penalty  method,  (1)  are  taken  to 
be  the  governing  equations  and  (2)  is  the  constraint.  Ex¬ 
pressing  (1)  in  the  variational  form  gives 


<$I  (u,  v;  6u,  5v)  =  JJ{  +  u~  +  v|^- 
3v 

lTt  +  u37  '  VS7J  '  o7  3y 


,  .  .  1  3P  , 

1  6u  +  F7  ST  6u 


+  +  +  V^7]  6v  +  7T-  IS  5v}  dx  dy' 


(10) 


where  represents  the  first  variation  of  the  functional  I 
with  respect  to  u  and  v  and  6u  and  6v  are  the  variations  of 
u  and  v  respectively.  Application  of  Green's  theorem,  a 
generalization  of  the  integration  by  parts  gives 


<5 1  (u,  v:  5u,  6v)  *  Jf{  [|£  +  u§£  +  v|H]  6u  +  [§J  +  ug  +  vg-  ) 


-  -L  6(^  +  |£)}  dx  dy  +  $  [t  <5u  +  t  iv]  ds. 
e0  3x  3y  x  y 


5v 
(ID 


P  P 

Here  t  *  —  n  ,  t  =  —  n  ,  and  (n  ,n  )  are  the  components 
x  p0  x  y  0o  y  x  y 

of  the  unit  normal  to  the  boundary  90.  The  integration  by 
parts  resulted  in  the  boundary  term  which  requires  that  eith¬ 
er  the  tractions  (t  and  t  )  or  the  variations  of  u  and  v  be 

x  y 

specified  on  90.  For  this  problem,  u  and  v  are  held  constant 
on  30  making  <5u  =  0  and  iv  *  0  on  the  boundary  30.  With  this 
specification  and  after  substituting  (2),  the  expression  (11) 
reduces  to 


61 (u, v; 6u, 6v) 


-jQF'ff  *  “i + 


+  l7£  +  US  +  vIt716v}  dx  dy- 


3v 


(12) 
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The  substitution  of  the  constraint  (2)  above  does  not  insure 


that  non-divergence  will  be  satisfied;  therefore,  it  is  added 
to  the  variational  form  (12)  as  a  weak  constraint.  This  is 
analogous  to  the  weak  constraint  method  described  by  Sasaki 
(l.)70a,b).  Adding  the  weak  constraint  to  (12)  gives 


<5l  (u,  v;  5u;  6v)  =  1 1|  +  u§£  +  v|^]  5u  +  [|£  +u§£  +vg]  6v 


+  X(l£  +  fy}  6u  +  Ty  Sv) }  dx  dy' 


where  X  is  the  penalty  parameter  (or  weight) . 

In  order  that  the  functional  I  attains  its  minimum  it 


is  necessary  that  6l  =  0.  It  should  be  noted  that  in  (13), 
the  pressure  has  been  eliminated  as  a  variable  so  that  the 
number  of  unknowns  has  been  reduced  to  two.  By  finding  the 
Euler-Lagrange  equations  for  (10)  and  (13)  and  then  equating 
the  like  terms,  one  can  find  the  following  relationship  for 
the  pressure  (Reddy,  1979b) . 


limit 


x  (^r  + 


-p. 


(14) 


This  is  based  on  the  fact  that  as  the  penalty  parameter  (X) 
approaches  infinity,  the  solution  (u^,  v^)  to  (13)  provides 
an  approximation  to  the  pressure  based  on  the  chosen  value 
for  X.  This  is  the  ’penalty’  resulting  from  the  elimination 
of  pressure  as  a  variable. 

For  the  penalty  finite  element  method,  the  velocity 
components  are  approximated  as  in  (6): 


(15) 


ue  =  £ue  n-  and  ve  =  Zve  N. . 
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substitution  of  (15)  into  (13)  gives  the  semidiscrete  Galer- 

kin  form 
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The  computational  form  for  (16)  is 


[Ce]  {A®}  +[(Ke]+X[pe]]  {A6}  =  {0}, 


(17) 
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For  computation,  matrices  [C  ]  and  [P  ]  need  be  computed  and 
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stored  for  each  element  only  one  time,  matrix  [K  ]  is  non¬ 
linear  and  must  be  recomputed  on  every  time  step  and/or  iter¬ 
ation. 


Finite  Element  Model  of  Shallow  Water  Equations 


The  models  described  above  provide  a  good  test  of  the 
ability  of  the  FEM  to  handle  the  nonlinear  effects,  which  are 
important  terms  in  the  governing  equations  of  the  atmosphere; 
however,  one  great  difficulty  in  operational  numerical  weath¬ 
er  prediction  is  the  occurence  of  fast-moving  waves  such  as 
gravity  waves.  Models  must  be  able  to  handle  these  waves 
without  their  energy  increasing  and  consequently  destroying 
the  forecast.  A  good  model  to  test  these  effects  is  the  well 
known  shallow  water  model.  The  shallow  water  equations  in  u 
and  v  are  the  same  as  (1)  except  that  the  Coriolis  term 
appears.  The  third  equation  which  governs  the  height  of  the 
free  fluid  surface  results  from  integration  of  the  contin¬ 
uity  equation  in  the  vertical  (Haltiner,  1971)  .  The  shallow 
water  equations  are  given  by 
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(18) 


Here  f  is  the  Coriolis  term,  h  is  the  height  of  the  free 
surface  of  the  fluid,  and  g  is  gravitational  acceleration. 
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The  semidiscrete  Galerkin  form  of  (18)  is 
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If  u,  v,  and  h  are  approximated  over  a  typical  element  by  u 
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then  (19)  becomes 
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The  only  matrix  in  (20)  which  does  not  change  is  [Me] .  All 
others  must  be  recomputed  for  each  element  on  each  iteration 


and/or  time  step. 


The  matrix  equations  (8),  (17),  and  (20)  are  a  result 
of  the  general  application  of  the  FEM.  At  this  point  a  gen¬ 
eral  computer  program  could  be  written  for  each  model.  Once 
the  general  program  is  developed,  consideration  must  be  given 
to  actual  space  and  time  discretization  which  includes  the 
development  of  the  element  interpolation  functions.  The  fin¬ 
al  steps  are  the  selection  of  the  actual  mesh  for  the  domain 
of  the  model  and  the  imposition  of  the  specific  initial  and 
boundary  conditions  required  for  the  solution.  The  discus¬ 
sion  which  follows  details  the  approach  taken  in  each  of  these 
areas  for  the  three  models  considered. 
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Numerical  Procedure 


Space  Discretization 

The  finite  element  method  involves  the  approximation 
of  functions  over  subdomains  (fie>  of  the  global  domain  (ft) . 
The  possible  shapes  of  these  subdomains  or  elements  are  lim¬ 
itless;  however,  unusual  or  completely  irregular  shapes  are 
difficult  to  manage  and  their  corresponding  interpolation 
functions  are  not  easily  obtainable.  The  most  widely  used 
element  shapes  in  two  dimensional  finite  element  analysis  are 
the  triangle  and  rectangle.  The  order  of  the  interpolation 
functions  used  determines  the  number  and  placement  of  nodal 
points  within  the  element  shape. 

In  this  study  five  different  elements  are  tested:  the 
three-node  linear  triangle,  the  four-node  bilinear  rectangle, 
the  six-node  quadratic  triangle,  the  eight-node  quadratic 
rectangle,  and  the  nine-node  quadratic  rectangle.  The  eight - 
and  nine-node  elements  both  use  quadratic  interpolation  but 
the  interpolation  functions  are  developed  in  different  ways. 
The  determination  of  the  interpolating  functions  is  discussed 
in  Appendix  A.  These  interpolating  functions  correspond  to 
the  term  shown  in  (6,  7,  8,  9,  15,  16,  19,  and  20)  so  that 
the  semi-discrete  form  of  the  equations  is  developed  indepen¬ 
dent  of  the  element  choice  and  is  then  general  for  any  ele¬ 
ment  application. 

For  the  problem  discussed  here,  the  element  size  is 
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held  constant  throughout  the  domain  of  the  problem.  This  has 
many  simplifying  effects  on  the  calculations  and  greatly  re¬ 
duces  computer  storage  requirements  since  element  matrices 
such  as  [Me]  in  (20)  need  be  calculated  for  only  one  element 
and  stored  rather  than  for  each  element.  In  a  model  with  com¬ 
plex  boundaries  or  where  grid  refinement  is  desired  this  sav¬ 
ings  is  not  possible. 

When  the  element  size  is  constant  it  is  possible  to 
reduce  the  FEM  to  a  finite  difference  form;  however,  the  re¬ 
sulting  finite  difference  form  will  be  valid  only  for  an  in¬ 
terior  grid  point  and  special  considerations  must  be  made  near 
the  boundaries.  Jesperson  (1974)  developed  a  finite  differ¬ 
ence  form  for  the  vorticity  -  stream  function  equation  (8) 
using  the  bilinear  rectangle  finite  element  approximations. 

The  resulting  finite  difference  equation  is  given  by 
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where  d  is  the  spacing  between  grid  points  and  the  subscripts 
indicate  the  grid  point.  As  noted  by  Jesperson,  the  right 
side  of  (21)  is  identical  to  the  well-known  Arakawa  Jacobian; 
however,  this  expression  for  the  Jacobian  was  developed  in  a 
completely  different  way  by  Arakawa  (1966) . 

In  this  study,  the  results  using  (21)  as  well  as  a 
simple  centered  space  finite  difference  scheme  are  compared 
to  the  results  using  the  five  different  finite  elements.  The 
left  side  of  (21)  will  be  discussed  in  a  later  section  deal¬ 
ing  with  consistent  and  lumped  mass. 

Time  Discretization 

Several  time  differencing  schemes  have  been  used  for 

% 

comparison  of  their  effects  on  the  solution.  The  Crank - 
Nicholson  scheme  was  chosen  as  the  basis  of  comparison  for  the 
different  finite  element  discretizations  because  it  has  neu¬ 
tral  stability,  no  computational  mode,  and  therefore  should 
have  no  effect  on  the  phase  speed  or  the  amplitude  of  the 
solution  (Mesinger  and  Arakawa,  1976) .  The  major  drawback  to 
this  scheme  is  that  it  is  fully  implicit  and  requires  itera¬ 
tion  at  each  time  step.  The  Crank-Nicholson  method  is  includ¬ 
ed  in  a  family  of  schemes  known  as  the  0 -family  of  approxima¬ 
tions  (Zienkiewicz ,  1977)  which  may  be  expressed  by 

(g}n+1  =  {q)n  +  9 At  {q}n+1  +  (1-6)  At{q}n  ,  (22) 
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where  the  superscript  indicates  the  time  (t=nAt) ,  (ql  is  the 
unknown  function,  At  is  the  time  step,  and  6  determines  the 
specific  scheme.  For  example,  some  of  the  possible  schemes 
are, 

9=0,  Forward  (Euler) , 

9  =  1/2,  Crank-Nicholson, 

9  =  2/3,  Galerkin,  and 
9=1,  Backward. 

Of  these  four  schemes,  all  were  used  except  the  Forward  dif¬ 
ferencing  scheme  which  was  eliminated  because  linear  stabil¬ 
ity  analysis  indicates  that  it  is  always  unstable.  All  of  the 
remaining  schemes  are  fully  implicit.  The  Galerkin  and  Back¬ 
ward  methods  both  damp  the  amplitude  of  the  solution  with  the 
greater  damping  occurring  in  the  Backward  scheme.  In  addition, 
two  explicit  schemes  will  be  used:  the  well-known  Leap-frog 
or  centered  scheme  which  is  given  by 

(q}n+1  =  (q}n_1  +  2At  {q}n  ,  (23) 

and  the  Matsuno  scheme  which  is  a  two-step  method 
given  by 

r  ,n+l*  {  ,n  .  ^  f',n 

tq}  =  {q}  +  At  { q } 

(24) 

(q>n+1  =  (q}n  *  At  tq>"+1*. 

the  Matsuno  scheme  results  in  damping  of  the  solution.  The 
Leap-frog  scheme  suffers  from  the  presence  of  the  computation¬ 
al  mode  which  restricts  the  choice  of  time  step. 
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This  restriction  may  be  determined  by  applying  the 
von  Neumann  stability  analysis  to  a  linearized  form  of  the 


advective  equation  (Mesinger  and  Arakawa,  1976).  For  a  cen¬ 
tered  finite  difference  scheme  applied  to  the  linear  version 
of  (5)  the  von  Neumann  technique  gives  the  well  known  CFL 
criteria  for  stability: 

n  Sp.  <  i,  (25) 


where  c  is  the  phase  speed  of  the  wave  and  the  factor  /2  is 
present  for  the  two-dimensional  problem.  When  the  FEM  is 
applied  to  the  linear  version  of  (5)  and  the  four-node  bi¬ 
linear  element  is  used  with  Leap-frog  time  differencing  the 
time  step  restriction  is 


(26) 


Thus,  the  FEM  imposes  a  penalty  of  / 3  on  the  time  step  (Cullen, 
1974).  The  time  step  penalty  results  because  of  the  coupling 
of  the  time  derivative  through  the  mass  matrix.  Some  detailed 
stability  analyses  are  given  in  Appendix  B. 

This  coupling  of  the  time  derivative  shown  by  the  left 
side  of  (21)  is  known  as  consistent  mass.  Physically,  (21) 
states  that  the  advection  which  is  calculated  by  the  right  side 
of  the  equation  for  the  node  (i,j)  affects  the  time  derivative 
not  only  at  that  node  but  also  at  the  surrounding  nodes.  In 
contrast,  when  lumped  mass  is  used  the  left  side  of  (21)  de¬ 
generates  to  a  single  term  {?.  .).  The  resulting  equation 
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is  the  traditional  approach  taken  in  finite  difference  methods 
and  as  shown  here  and  by  Sasaki  and  Chang  (1979) ,  the  results 
using  (21)  are  much  more  accurate  than  the  traditional  finite 
difference  approach.  As  discussed  by  Zienkiewicz  (1977) 
there  are  several  approaches  to  lumping  the  mass  in  the  FEM; 
however,  the  simplest  method  is  to  make  the  matrix  [Me]  in 
(8),  (17),  and  (20)  diagonal.  The  latter  method  is  used  in 
this  study.  The  choice  of  lumped  or  consistent  mass  greatly 
affects  the  computational  details  of  the  problem  being  con¬ 
sidered.  Computational  details  will  be  discussed  in  the 
following  section. 


Computational  Details 

The  vorticity  equation  (8)  combined  with  Crank-Nich- 
olson  time  differencing  may  be  written  as  either 


or  [M]  {A?}=  {F}  ,  (28) 

where  {f}  =  -  ^  [J]n+1  U }n+1  +  [[M]  -  ^  [J]n]  U)n  , 

{C}n+1  =  {£}n  +  {A?}. 

From  an  algebraic  point  of  view  (27)  and  (28)  are  equivalent 
equations;  however,  computationally  they  are  quite  different. 
When  solving  (27),  the  right  side  of  the  equation  is  known. 

The  matrix  [J]n+^ 


is  not  known  since  it  depends  on  the  stream 


function  at  the  new  time.  The  procedure  for  solving  (27)  is 
to  calculate  the  right  side,  then  calculate  [J]n+1  based  on 
a  guess,  then  solve  for  The  stream  function  equation 

(9)  is  then  solved  for  {<Jj}n+^.  These  new  values  of  iJj  serve 
as  the  guess  in  the  following  iteration.  The  iterative  pro¬ 
cess  is  repeated  until  the  solution  converges.  The  coeffi¬ 
cient  matrix  on  the  left  is  recomputed  for  each  element  on 
each  iteration  and  assembled  into  the  global  matrix  which  is 
banded  and  unsymmetric.  The  result  is  a  set  of  linear  alge¬ 
braic  equations  which  must  be  solved  for  {£}n+^.  The  latter 
task  is  extremely  time  consuming  for  the  computer,  especially 
when  the  number  of  nodes  is  large. 

_  i  I 

The  alternative  is  to  use  (28).  The  matrix  [j] 
must  still  be  computed  but  it  is  multiplied  by  a  guess  (values 
from  the  previous  iteration)  for  {^}n+1  and  placed  into  the 
force  vector.  The  coefficient  matrix  [M]  is  banded  and  sym¬ 
metric.  It  must  be  calculated  only  once  and  can  be  decom¬ 
posed  once  using  Cholesky  decomposition  (Carnahan,  et  al . , 
1969) .  Then  on  each  iteration,  only  forward  and  backward  sub¬ 
stitution  are  required  to  calculate  the  solutions.  Even 
though  the  number  of  iterations  required  for  convergence  may 
be  greater  for  (28)  than  for  (27),  the  overall  time  savings 
can  be  substantial. 

When  using  an  explicit  time  differencing  method  the 
resulting  form  is  similar  to  (28)  and  the  set  of  algebraic 
equations  can  be  solved  as  in  (28)  with  no  iteration.  The 


latter  point  makes  the  Leap-frog  differencing  scheme  very 
attractive  computationally.  If  lumped  mass  is  used  the 
matrix  [M]  becomes  diagonal  and  the  solution  of  (28)  is 
relatively  simple  and  fast.  The  combination  of  lumped  mass 
and  (27)  results  in  little,  if  any,  savings.  Specific  re¬ 
sults  from  each  method  will  be  discussed  in  Chapter  IV. 


The  Model  Problem 

For  the  vorticity-stream  function  model  the  domain 
chosen  is  a  channel,  3800  km  in  length  and  width. 

The  initial  conditions  for  the  stream  function  are 

given  by 

ip*  (x,  y)  -  -V0y  -  exp [-a  (x2+y2 )  ] ,  (29) 


where  U0,  ,  and  a  are  constants  and  the  origin  is  at  the 

center  of  the  domain.  This  model  problem  is  identical  to  the 
one  used  by  Sasaki  and  Reddy  (1979)  and  Sasaki  and  Chang  (1979) 
and  is  used  here  to  facilitate  comparison  of  results.  The 
initial  conditions  for  vorticity  were  derived  by  appropriate 
differentiation  of  (29) .  These  initial  conditions  are  shown 
in  Fig.  1.  The  boundary  conditions  are 
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east-west  boundaries 
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where  the  channel  length  is  L,  w  is  the  width,  and  and  C2 
are  constants.  For  these  initial  and  boundary  conditions, 
the  model  equations  (5)  have  a  unique  analytical  solution. 

y(x,y,t)  =  <l>*  (x-U0t ,  y)  .  (30) 

The  solution  indicates  that  the  circular  vortex  propagates 
eastward  until  it  reappears  on  the  west  edge  and  returns  to 
its  initial  position  after  a  period  of  L/U0.  For  the  con¬ 
stants  chosen  here,  that  period  is  exactly  120  hours. 

For  the  penalty  method,  the  initial  conditions  are 
obtained  by  using  the  definition  (3)  and  the  initial  condi¬ 
tions  for  stream  function  (29) .  The  resulting  initial  con¬ 
ditions  for  u  and  v  describe  the  same  circular  vortex.  These 
initial  conditions  are  shown  in  Fig.  2.  As  long  as  the  proper 
value  for  the  penalty  parameter  (X)  is  chosen,  the  behavior 
of  the  vortex  should  be  similar  to  the  behavior  experienced 
in  the  vorticity  -  stream  function  model.  The  boundary  con¬ 
ditions  used  are 

u  (0 ,  y)  =  u  (L, y)  ,  *1 

\  east-west  boundaries 

v (0 , y)  =  v(L, y) ,  J 

u(x,o)  =  u(x,W)  =  U0, 

l  north-south  boundaries 

v(x,o)  =  v(x,W)  =  0.  J 

Here  the  u  and  v  components  must  both  be  specified  because  of 
the  boundary  term  which  appears  in  (11) . 

In  the  shallow  water  equation  model,  the  domain  and 
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initial  conditions  chosen  are  the  same  as  used  by  Grammelt- 
vedt  (1969)  and  later  by  Cullen  (1974a) .  The  domain  is  a 
channel  6000  km  in  length  and  4400  km  wide.  The  initial 
conditions  are 


h  ( x ,  y ) 


where 


Ho+Hl  *  tanh  +H2  •  sech^iiXZYoi  . 

(30) 

[0.8  sin  (~)  +  0.5  sin  (^L)  ]  , 

Li  Lt 


H0  =  2000  m 
H,  =  -22  0  m 
H2  =  133  m, 
g  =  10  ms'2 


L  =  6000  km 
W  =  4  4  00  km 
Y0=  W/2, 


The  initial  conditions  for  u  and  v  are  determined  using  the 
FEM  with  the  geostrophic  relation, 

,,  _  g  3h  _  g  9h 

u  '  ¥  3y  and  v  -  t  7Z  ’ 


The  'Beta-plane'  approximation  is  made  where-in  the  Coriolis 
parameter  (f)  is  calculated  from 


f  -  f.  +  By, 

where  fo=1.0*10  ^  and  6  =  1.5*10  ^m  ^ .  These  initial 

conditions  are  shown  in  Fig.  3.  Physically,  these  initial 
conditions  describe  a  west  to  east  jet-stream  which  has  north- 
south  disturbances  along  its  axis  (Grammeltvedt ,  1969) .  The 
initial  conditions  are  shown  in  Fig.  3.  The  boundary  condi- 
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tions  specified  for  this  model  are 
u (0, y)  =  u(L,y) , 

v(0,y)  =  v(Li,y),  \  east-west  boundary 

h(0,y)  =  h (L, y) ,  1 

v(x,0)  =  v (x, W)  =  0.  north-south  boundary 

It  is  interesting  to  note  that  there  is  a  subtle  difference 
between  the  north-south  boundary  conditions  specified  here 
and  those  required  in  finite  difference  models.  Grammeltvedt 
(1969)  was  forced  to  specify  boundary  conditions  for  all  var¬ 
iables  in  his  models;  this  is  an  overspecification  in  a  con¬ 
tinuous  model  and  sometimes  called  numerical  boundary  condi¬ 
tions.  This  additional  specification  at  the  boundaries  can 

tt 

adversely  effect  the  solution  (Sundstrom,  1973) . 

The  Grid  (Finite  Element  Mesh) 

Both  the  vorticity-stream  function  and  penalty  models 
are  solved  on  an  evenly-spaced  grid  consisting  of  11  x  11 
nodal  points  with  380  km  spacing.  For  the  shallow  water  equa¬ 
tions,  a  grid  of  15  x  21  nodal  points  is  employed  with  300  km 
spacing  in  the  east-west  direction  and  approximately  315  km 
spacing  in  the  north-south  direction.  In  order  to  impose  per¬ 
iodic  boundary  conditions  in  the  FEM,  the  nodal  values  on  the 
eastern  boundary  are  considered  to  be  identical  to  those  on 
the  western  boundary  but  their  locations  are  preserved. 

As  stated  earlier,  linear  and  quadratic  rectangular 
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as  well  as  triangular  elements  are  tested  for  their  relative 
accuracy.  Regardless  of  the  element  type  (i.e.,  rectangular 
with  four-or  nine-nodes,  or  triangular  with  three-or  six- 
nodes)  ,  the  number  of  elements  chosen  is  such  that  the  total 
number  of  nodes  remains  the  same.  That  is,  one  nine-node 
quadratic  element  replaces  four  bilinear  four-node  elements; 
or,  two  three-node  linear  triangular  elements  replace  one 
four-node  linear  element,  etc.  However,  when  an  eight-node 
quadratic  (serendipity)  element  is  used  (the  center  node  is 
missing  when  compared  to  the  nine-node  element) ,  the  total 
number  of  nodes  is  reduced.  In  the  case  of  the  vorticity- 
stream  function  model,  the  total  number  of  nodes  is  reduced 
from  121  to  96.  Fig.  4  shows  the  relative  size  and  arrange¬ 
ment  of  the  elements  used.  Both  the  three-and  six-node  tri¬ 
angles  shown  in  Fig.  5  have  been  tested  with  the  vorticity- 
stream  function  model.  These  results  are  presented  in  Chap¬ 
ter  IV. 

Accuracy  Indicators 

The  results  from  the  vorticity-stream  function  and 
penalty  models  using  the  different  time  and  space  discretiza¬ 
tions  discussed  are  compared  for  lumped  versus  consistent 
mass,  position  of  the  nonlinear  terms,  convergence  rate,  pro¬ 
gram  size,  computer  time,  and  accuracy  of  the  solution.  Other 
than  accuracy,  all  the  comparisons  are  straight  forward.  The 
root-mean-sqaure  error  (RMSE) ,  given  by 
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Fig.  5.  Arrangement  of  triangular  elements  within 
the  mesh 


RMSE 


Z (Fi_Si) ** 

N 

where  N  is  the  number  of  grid  points,  are  the  forecast 
values,  and  are  the  solution  values,  is  used  as  the  basis 
of  comparison  for  accuracy.  This  measure  can  be  misleading 
if  used  alone.  Here  the  RMSE  is  combined  with  visual  examin¬ 
ation  of  the  analytical  solution  and  120  hour  forecast  fields 
of  the  vorticity  and  stream  function. 

Comparisons  of  the  results  from  the  shallow  water 
equation  model  are  more  complicated  since  the  equations  are 
non-linear  and  there  is  no  analytical  solution.  However, 
there  are  some  facts  about  the  model  which  can  be  checked 
and  compared. 

The  absolute  vorticity,  and  the  total  energy  are  con¬ 
served  in  the  FEM  model .  Another  parameter  conserved  in  the 
FEM  model  is  the  mean  wave  number  (Cullen,  1974b) .  Because 
of  this,  the  allocation  of  energy  among  wave  numbers  should 
not  change  significantly  during  the  forecast.  Therefore,  a 
two  dimensional  Fourier  analysis  has  been  performed  on  the 
initial  height  fields  as  well  as  the  48~hour  forecast  height 
fields.  These  analyses  are  based  on  a  discussion  by  Goodman 
(1968)  and  calculated  using  an  algorithm  developed  by  Duchon* 
Fig.  6  shows  the  two-dimensional  Fourier  analysis  of  the  in¬ 
itial  height  field.  This  three  dimensional  plot  shows  the 

■''Associate  Professor  C.E.  Duchon,  School  of  Meteor¬ 
ology,  University  of  Oklahoma  provided  a  copy  of  the  algorithm 
for  this  work. 
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shallow  water  equations 


relative  amplitudes  for  the  combination  of  wave  numbers  in 
the  x  (Along  the  channel)  direction  and  y  (across  the  chan¬ 
nel)  direction.  Basically,  the  height  field  slopes  from 
south  to  north  requiring  all  wave  numbers  across  the  channel 
to  be  represented.  Only  wave  numbers  1  and  3  are  represented 
along  the  channel.  The  amplitudes  are  plotted  on  a  logrith- 
mic  scale  so  that  the  smaller  amplitudes  will  be  visible. 


CHAPTER  IV 

RESULTS 

As  described  earlier,  the  three  models  used  were 

1.  Vorticity  -  Stream  Function, 

2.  Penalty  Method,  and 

3.  Shallow  Water  Equations. 

The  tests  performed  using  these  models  answer  the  questions 
relating  to  element  choice,  selection  of  time  differencing 
scheme,  and  use  of  lumped  or  consistent  mass.  Here  the  test 
results  are  presented  for  each  of  the  models.  Following  in 
Chapter  V  will  be  a  discussion  of  the  conclusions  which  may 
be  drawn  from  these  test  results  taken  as  a  whole. 

Vorticity  -  Stream  Function  Model 

Table  1  shows  the  results  for  the  comparison  of  the 
five  elements  tested  using  Crank-Nicholson  time  differencing 
and  consistent  mass.  The  elements  are  compared  as  to  position 
of  the  nonlinear  terms,  convergence  rate,  program  size,  com¬ 
puter  time,  and  accuracy.  The  use  of  (27)  with  the  nonlinear 
terms  on  the  left  side  cr  (28)  where  they  are  on  the  right 
makes  no  difference  in  the  accuracy.  For  each  element,  this 
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Table  1.  Comparison  of  5  Finite  Elements  using  Crank-Nicholson  Time-Difference 


result  is  the  same.  However,  the  computer  time  required  when 
the  non-linear  terms  are  on  the  left  is  nearly  three  times 
as  great.  In  addition,  there  is  a  savings  of  computer  stor¬ 
age  when  the  nonlinear  terms  are  on  the  right  side.  The  only 
drawback  r.o  moving  the  nonlinear  terms  to  the  right  side  is 
the  slight  increase  in  the  average  number  of  iterations  re¬ 
quired  for  convergence  on  each  time  step  in  the  quadratic 
elements . 

Fig.  7  shows  the  fields  after  120  hours  of  forecast 
for  the  four-and  three-node  element  using  consistent  mass. 

The  analytic  solution  is  presented  for  comparison.  Fig.  8 
gives  the  same  information  for  the  six-and  nine-node  elements. 
All  of  the  elements  give  reasonable  solutions;  however,  the 
four-node  result  is  clearly  the  best.  The  phase  speed  for 
the  three-and  four-node  elements  is  very  accurate  but  the 
amplitude  is  reduced  slightly  more  in  the  three-node  than  the 
four-node  solution.  The  nine-node  solution  has  a  phase  speed 
which  is  too  small  while  the  six-node  solution  has  phase 
speed  which  is  slightly  high.  Both  methods  result  in  some 
distortion  of  the  field  with  the  presence  of  noise  evident  in 
the  six-node  vorticity  field.  The  results  from  the  eight-node 
element  (not  pictured)  are  very  poor.  The  visual  comparisons 
are  supported  by  the  RMSE  given  in  Table  1.  Not  only  is  the 
RMSE  the  lowest  for  the  four-node  element  but  the  program  size 
and  computer  time  are  also  smaller. 


48 


jinrun  rune!  ia«  «■/ r, 


X 

ta 

fa 


<U 

•O 

0 

Z 

I 

ro 


C 

o 

•H 

P 

3 


O 

W 


« 

c 

< 


X 

w 

fa 


P 

c 

(0 


o 

•P 

P 


0) 

T3 

O 

z 

( 


c 

o 


<P 


<u 
<u  CJ 


TJ 

o 

c 

I 


e 

id 
<U 
p 
P 
W 
I 

>1 
P 
•H 

o 

e  >h 

<U  P 


Tf 

C 

id 


n  u-i 
p 


r)  -h  <u 


c  -a  x: 


m 

i 

p- 


P  P  P 


P  JS 


U] 
p 
3 
O 
O  X 


id 


o  o 

a-H  (n 

e  z  p 

c  I 

U  M  p 
c  <u 
aj  p 
.  p  <p 

Is"  U  id 


•  JS  to 
tji  p  to 
•p  -p  id 
hie 


model . 


Fig.  8.  Comparison  of  6-and  9-node  finite  elements  with  Crank-Nicholson 
time  difference  and  consistent  mass  after  120  hours  for  the  vorticity- 
stream  function  model. 


As  discussed  in  Chapter  III  there  are  several  possible 
configurations  of  both  the  three-and  six-node  elements  when 
the  nodal  spacing  is  held  constant.  The  triangular  arrange¬ 
ment  chosen  does  effect  the  solution.  Table  2  gives  the  RMSE 
for  each  of  the  five  configurations  (shown  in  Fig.  5)  tested 
for  both  the  three-and  six-node  triangular  elements.  For 
both  elements,  the  case  in  which  the  triangular  slant  alter¬ 
nates  in  the  direction  of  flow  gives  the  lowest  RMSE .  A  vis¬ 
ual  examination  of  the  120-hour  forecast  for  the  three-node 
element  (Fig.  9a-d)  reveals  little  difference  in  the  fields. 
The  vorticity  pattern  in  Fig.  9a, b  is  slightly  distorted  in 
the  direction  opposite  to  the  constant  diagonal  slant  but 
this  distortion  is  not  evident  in  the  six-node  element  (Fig. 
lOa-d) .  In  the  six-node  element  the  case  which  has  constant 
diagonal  slant  from  lower  left  to  upper  right  (Fig.  10a)  has 
the  most  accurate  phase  speed  and  its  RMSE  is  not  significant¬ 
ly  different  from  Fig.  lOd  which  has  alternating  slant  in  the 
direction  of  flow.  It  is  evident  that  the  linear  triangle 
handles  the  advection  much  more  accurately  than  the  quadratic 
regardless  of  the  triangle  orientation.  The  diagonal  slant 
is  examined  further  for  the  three-node  element  using  the 
shallow  water  equations. 

Table  1  is  based  on  the  use  of  the  consistent  mass 
matrix  with  Crank-Nicholson  time  differencing.  In  contrast, 
Table  3  contains  the  same  comparison  of  elements  for  the  lump¬ 
ed  or  diagonalized  mass  matrix.  The  RMSE  is  considerably 


Table  3.  Comparison  of  5  Finite  Elements  using  Crank -Nicholson  Time  Difference  with 
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higher  when  the  mass  is  lumped,  regardless  of  the  element 
used.  Fig.  11  graphically  depicts  this  increase  in  RMSE  for 
the  three-and  four-node  elements.  The  phase  speed  is  reason¬ 
ably  accurate  but  the  vorticity  center  moves  southward  with 
considerable  loss  of  amplitude.  As  shown  in  Fig.  12,  the  six 
and  nine-node  elements  which  use  quadratic  interpolation  pro¬ 
vide  the  better  results  when  the  mass  is  lumped.  It  is  evi¬ 
dent  in  the  nine-node  element  that  most  of  the  RMSE  is  due  to 
inaccuracy  in  the  phase  speed. 

In  summary,  the  four-node  bilinear  element  has  result 
ed  in  the  most  accurate  120-hour  forecast.  The  phase  speed 
and  amplitude  are  close  to  the  analytical  solution  and  the 
RMSE  is  the  lowest.  This  element  is  used  to  compare  five 
time  differencing  schemes  using  both  consistent  and  lumped 
mass.  These  results  appear  in  Table  4  (includes  comparative 
information  on  four-node  element  repeated  from  Table  1  and  3) 
After  120  hours  the  RMSE  for  Leap-frog  differencing  is  the 
same  as  for  Crank-Nicholson.  Further,  it  is  difficult  to 
detect  any  difference  between  Fig.  13  which  shows  the  Leap¬ 
frog  result  and  Fig.  7  which  shows  the  Crank-Nicholson  result 
The  other  difference  schemes  all  give  reasonably  low  RMSE. 

As  shown  in  Figs.  13  and  14,  the  Matsuno,  Galerkin,  and  Back¬ 
ward  methods  are  all  accurate  in  phase  speed  but  damping  of 
the  amplitude  accounts  for  the  increased  RMSE.  Figs.  15  and 
16  show  120-hour  forecasts  for  the  same  time  differencing 
schemes  using  lumped  mass.  These  lumped  mass  results  do  not 
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16.  Comparison  of  Galerkin  and  Backward  time  differencing  with 
iped  mass  after  120  hours  using  the  4-node  finite  element  for  the 
ticity-stream  function  model. 


differ  greatly  from  the  previous  discussion  nor  from  the  re¬ 
sults  of  Sasaki,  and  Chang,  (1979)  . 

Table  5  gives  the  results  from  some  finite  difference 
schemes  for  comparison.  The  Arakawa  difference  scheme  with 
consistent  mass  is  given  in  (21)  .  As  stated  previously,  this 
scheme  is  equivalent  to  the  four-node  finite  element  and  in¬ 
deed,  the  results  are  the  same.  The  only  difference  in  the 
two  schemes  is  that  no  north  and  south  boundary  conditions 
need  be  specified  for  vorticity  in  the  FEM  but  they  are  nec¬ 
essary  to  carry  out  the  finite  differencing.  In  this  pro¬ 
blem,  there  was  no  noticeable  difference  in  the  result;  how¬ 
ever,  this  may  not  always  be  the  case  (Sundstrom,  1973)  .  It 
should  be  noted  that  for  a  particular  problem,  when  the  el¬ 
ements  are  a  regular  shape  and  of  equal  size  everywhere,  the 
FEM  can  be  simplified  to  a  quasi-FD  model  with  corresponding 
savings  in  computer  time  and  storage.  Here  the  time  require¬ 
ment  decreased  three-fold  and  the  storage  requirement  was 
half  (Table  5) .  Figs.  17  and  18  show  the  result  for  both 
consistent  and  lumped  mass  for  Arakawa  and  centered  finite 
difference  methods.  The  traditional  finite  difference  prac¬ 
titioner  would  achieve  results  similar  to  those  with  lumped 
mass.  With  consistent  mass,  both  the  Arakawa  and  centered 
difference  methods  give  good  results.  As  shown  in  Table  5, 
the  use  of  consistent  mass  results  in  approximately  a  40 
percent  increase  in  computer  time  over  the  traditional  ap¬ 
proach. 
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Finally,  Table  6  shows  the  results  when  the  three 
best  schemes  are  intergrated  for  a  longer  time  using  differ¬ 
ent  time  steps.  As  discussed  earlier,  1  hour  is  the  upper 
limit  for  the  Leap-frog  integration.  Other  than  the  Crank- 
Nicholson  method  with  2-hour  time  step,  the  result  is  not 
significantly  different  to  240  hours.  By  360  hours  the 
Leap-frog  scheme  is  affected  by  its  computational  mode  but 
the  Crank-Nicholson  method  allows  extremely  long  integrations. 
After  1080  hours,  the  Crank-Nicholson  method  finally  gives  a 
fairly  high  RMSE;  however,  this  is  due  mainly  to  phase  speed 
error.  The  shape  of  the  fields  is  quite  good. 

Penalty- Method 

The  testing  conducted  using  the  penalty  method  is 
much  less  extensive  than  for  the  vorticity-stream  function 
model.  Table  7  gives  the  results  from  comparison  of  elements 
using  Crank-Nicholson  time  differencing.  The  lowest  RMSE 
occurs  for  the  four-node  bilinear  element.  It  is  followed  by 
the  nine-,  six-,  and  eight-node  elements.  The  order  is  ex¬ 
actly  the  same  as  in  the  previous  discussion.  Fig.  19  shows 
the  initial  and  120-hour  forecast  for  the  u  and  v  wind  com¬ 
ponents  as  well  as  the  stream  function  and  vorticity  calcula¬ 
ted  from  the  forecast  winds.  The  slight  loss  of  symmetry  in 
the  forecast  u  and  v  fields  indicates  that  the  divergence  free 
condition  is  not  exactly  satisfied.  In  addition,  the  model 
has  given  a  reduction  in  the  gradients  of  u  and  v  which  re- 
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suits  in  a  decrease  in  the  strength  of  the  vorticity  center. 

The  phase  speed  is  as  accurate  as  it  was  in  the  vorticity- 
stream  function  model  but  there  is  a  slight  northward  shift 
in  the  vorticity  pattern  which  was  not  present  in  the  vortic- 
ity-stream  function  model  when  the  four-node  element  was 
used  (Fig.  7).  From  Table  6,  it  is  evident  that  overall,  the 
penalty  method  is  not  as  accurate  for  this  advection-domin- 
ated  flow  situation  as  the  vorticity-stream  function  model. 

As  discussed  in  Chapter  III,  approximate  values  for 
the  pressure  may  be  calculated  using  (14).  The  calculated 
pressure  depends  not  only  on  the  accurate  calculation  of 
divergence  but  also  on  the  proper  selection  of  the  penalty 
parameter  X.  This  parameter  can  be  found  through  trial  and 
error.  As  discussed  by  Reddy  (1979b) ,  too  large  a  value  of 
X  results  in  the  degeneration  of  the  governing  equations  into 
only  the  continuity  equation;  too  small  a  value  means  that 
the  continuity  equation  will  not  be  satisfied.  Even  when  the 
best  value  for  X  is  found,  the  continuity  equation  is  not 
exactly  satisfied.  The  resulting  pressure  values  are,  at  best, 
an  approximation. 

The  benefit  from  the  use  of  the  penalty  method  is  the 
reduction  of  the  number  of  unknowns,  thereby  reducing  the 
number  of  equations  to  two.  The  form  of  the  equations  given 
by  (17)  is  the  traditional  approach  taken  in  the  application 
of  the  FEM.  The  computational  times  given  in  Table  6  using 
that  approach  are  very  long  when  one  considers  that  the  grid 


only  contains  11  x  11  points.  Several  attempts  have  been 
made  to  improve  the  computation  time,  such  as  moving  the 
nonlinear  terms  and  the  penalty  term  to  the  force  vector,  and 
separating  the  u  and  v  equations  in  order  to  reduce  the  order 
of  the  system  of  linear  equations  which  must  be  solved. 

Every  such  test  resulted  in  reduced  accuracy. 

It  is  because  of  the  inaccuracies  in  the  pressure  and 
the  expense  of  computation  that  the  penalty  method  is  not 
tested  to  a  greater  extent.  It  is  interesting  however,  that 
the  different  approach  taken  in  the  penalty  method  results  in 
the  same  relative  accuracy  of  the  elements  tested. 

Shallow  Water  Equations 

Eqs . (13)  used  for  this  model  present  a  challenge  for 
any  type  of  finite  discretization  because  they  allow  fast 
moving  gravity  waves.  The  gravity  waves  are  initially  of 
much  lower  amplitude  than  the  long  waves  but  they  can  grow 
and  after  several  time  steps  they  can  obscure  the  features 
of  interest  in  the  solution.  Much  effort  has  been  directed 
toward  the  control  of  gravity  waves  (e.g.  Kwizak  and  Robert, 
1971);  however,  since  the  intent  of  this  study  is  to  find  the 
best  finite  element  discretization,  the  gravity  waves  are 
allowed  to  develop  without  any  control  so  that  the  discreti¬ 
zation  may  be  found  which  most  efficiently  discourages  their 
growth. 

All  of  the  elements  previously  discussed  were  tested 


73 


except  the  eight-node  serendipity  element  which  was  elimin¬ 
ated  because  of  its  poor  performance  with  the  vorticity- 
stream  function  and  penalty  models.  The  comparison  of  run 
times  and  program  size  for  the  models  tested  is  given  in 
Table  8.  The  three  equations  were  solved  one  at  a  time  in 
each  time  step  and/or  iteration  (depending  on  the  time  dif¬ 
ference  scheme) .  As  occurred  in  the  vorticity-stream  func¬ 
tion  problem,  the  quadratic  (6  and  9  nodes)  elements  require 
much  more  computer  time  and  storage  than  the  three-and  four- 
node  elements.  The  difference  in  these  statistics  between 
the  three-and  four-node  elements  is  insignificant.  Other 
than  these  two  comparisons  there  are  no  other  quantitative 
comparisons;  however,  there  are  qualitative  differences  in 
the  results  produced  by  the  elements. 

Fig. ' 20  shows  the  48-hour  forecast  (144  time  steps  of 
20  minutes  each)  of  the  height  fields,  from  the  six-and  nine- 
node  elements  using  Crank-Nicholson  time  differencing.  Also 
shown  are  the  corresponding  variance  plots  from  the  Fourier 
analyses  of  those  fields.  Initially  (Fig.  3)  there  are  three 
equally  spaced  troughs  in  the  channel .  The  western  trough 
is  shallow  to  the  north  and  deeper  to  the  south.  The  cen- 
teral  and  eastern  troughs  are  sharper  and  deeper  to  the  north. 
After  48  hours,  the  troughs  move  approximately  one-third  the 
channel  length  to  the  east.  As  shown,  both  the  six-and  nine- 
node  elements  produce  roughly  the  same  trough  positions. 
Movement  is  slightly  faster  for  the  nine-node  case.  It  can 
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readily  be  seen  that  the  six-node  element  has  resulted  in  the 
generation  of  many  short  waves  at  almost  every  possible  wave 
number,  particularly  at  the  highest  wave  numbers.  In  con¬ 
trast,  the  nine-node  element  has  not  led  to  the  production 
of  any  short  waves.  Both  elements  have  produced  wave  number 
two  and  wave  number  five.  This  may  be  a  natural  consequence 
of  this  method  since  it  conserves  the  average  wave  number. 

Fig.  21  shows  the  48  hour  height  field  and  its  Four¬ 
ier  analysis  for  the  four-node  rectangular  element  using 
Crank-Nicholson  time  differencing.  The  position  of  the 
troughs  is  very  close  to  the  nine-node  result  but  the  western 
trough  is  not  as  deep  as  the  nine-node  case.  The  Fourier 
analyses  of  the  nine-and  four-node  elements  are  different  in 
two  important  respects.  For  the  four-node'  element,  the  var¬ 
iance  at  wave  number  five  along  the  channel  is  approximately 
60  percent  of  the  nine-node  result.  The  same  difference  is 
true  for  the  wave  associated  with  wave  numbers  two  and  four 
in  the  x  and  y  directions  respectively.  The  slightly  better 
result  for  the  four-node  element  is  consistent  with  the  vor- 
ti^ity-stream  function  and  penalty  models  and  is  significant 
when  one  considers  the  difference  in  computational  time  be¬ 
tween  the  two  elements.  Some  comparable  results  were  also 
achieved  with  the  three-node  element. 

The  three-node  element  is  used  to  show  the  effect  of 
diagonal  slant  of  the  triangles.  Figs.  22  and  23  show  the 
height  fields  after  48  hours  using  Crank-Nicholson  time 
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differencing  for  triangle  configurations  1-4  (shown  in  Fig. 

5)  and  the  corresponding  Fourier  analyses.  The  trough  pos¬ 
itions  are  all  the  same;  however,  their  character  differs 
depending  on  the  triangle  slant.  For  slants  1  and  2,  the 
features  tilt  slightly  in  the  opposite  direction  of  the  con¬ 
stant  slant.  This  is  analogous  to  a  distortion  found  in  the 
vorticity  pattern  in  Fig.  9.  The  height  fields  for  slants  1 
and  2  definitely  have  a  smoother  appearance  than  for  slants 
3  and  4.  This  is  evident  in  the  Fourier  analyses  of  the 
fields.  The  alternating  diagonal  slant  appears  to  encourage 
early  development  of  gravity  waves.  Of  these  four,  the  worst 
case  results  when  the  slant  alternates  across  the  direction 
of  flow.  This  is  consistent  with  the  vorticity-stream  func¬ 
tion  result  for  the  same  triangle  orientation.  The  constant 
diagonal  slant  cases  appear  to  give  the  better  results  for 
triangles. 

Compared  to  the  four-node  element  (Fig.  21)  the  result 
for  slant  2  is  almost  the  same.  For  slant  2,  there  is  higher 
variance  in  the  combination  of  wave  number  two  along  the  chan¬ 
nel  and  wave  numbers  four  and  five  across  the  channel  than 
there  is  in  the  four-node  case.  There  is  also  a  development 
in  wave  number  four  along  the  channel  which  is  not  present  in 
the  four-node  case.  Therefore,  it  appears  that  the  four-node 
element  has  slightly  better  performance  in  this  case. 

The  four-node  element  was  tested  with  four  other  time 
differencing  schemes  besides  Crank-Nicholson: 
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1)  Galerkin 


2)  Backward, 

3}  Matsuno,  and 
4)  Leap-frog. 

As  shown  in  Table  8,  all  of  these  schemes  were  run  with  the 
equations  formulated  so  that  all  terms  are  in  the  force  vec¬ 
tor  except  the  time  derivative.  Comparing  the  two  four-node 
Crank-Nicholson  cases  one  can  see  that  this  formulation  re¬ 
sults  in  a  five-fold  decrease  of  computer  time.  Because  no 
iteration  is  required,  the  explicit  schemes  are  much  faster. 

The  results  from  the  Galerkin,  Backward,  and  Matsuno 
schemes  are  all  similar.  All  of  these  schemes  cause  some 
deunping  of  the  waves.  The  48-hour  forecast  height  fields  for 
Matsuno  and  Galerkin  schemes  as  well  as  their  Fourier  ana¬ 
lyses  are  shown  in  Fig.  24.  Compared  to  the  Crank-Nicholson 
result  (Fig.  21)  both  of  these  schemes  give  a  smoother  fore¬ 
cast  field.  The  only  difference  in  the  Fourier  analyses  is 
the  lower  amplitude  for  wave  number  five  along  the  channel  for 
both  the  Galerkin  and  Matsuno  methods  compared  to  the  Crank- 
Nicholson  results. 

The  fastest  scheme  tested  is  the  Leap-frog  or  centered 
time  difference.  It  requires  only  one-fourth  the  computation 
time  of  Crank-Nicholson  and  as  shown  in  Fig.  25,  it  provides 
reasonable  results.  At  48  hours,  there  is  more  generation  of 
wave  number  five  along  the  channel  as  well  as  higher  cross 
channel  wave  numbers.  The  location  of  the  troughs  is  very 


Fig.  24.  48  hour  forecast  height  fields  and  Fourier  analyses  from 

shallow  water  equations  using  Galerkin  and  Matsuno  time  differencing. 


close  to  the  Crank -Nicholson  result. 

The  Leap-frog  and  Crank-Nicholson  schemes  are  used 
to  demonstrate  the  conservative  properties  of  the  finite 
element  method.  As  stated  earlier,  for  the  shallow  water 
equations,  both  the  total  energy  and  absolute  vorticity  are 
conserved.  Figs.  26  and  27  show  the  total  energy  and  vor¬ 
ticity  plots  for  both  schemes.  As  expected,  the  Leap-frog 
scheme  with  its  computational  mode  becomes  computationally 
unstable  earlier  than  Crank-Nicholson.  The  beginning  of  in¬ 
stability  is  evident  when  one  looks  at  changes  in  the  kinetic 
energy,  a  quantity  which  is  not  conserved.  For  the  Leap-frog 
scheme,  the  kinetic  energy  increases  less  than  three  percent 
for  the  first  thirty  days  and  then  rapidly  increases  near  the 
thirty-five  day  mark  (see  Fig.  28).  As  shown  in  Fig.  28,  the 
increase  in  kinetic  energy  is  much  slower  for  the  Crank- 
Nicholson  scheme. 
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ig.  28.  Kinetic  energy  (normalized)  versus  time  (days) 
or  Leap-frog  and  Crank-Nicholson  schemes  using  4 -node 
lement  with  shallow  water  equations. 


CHAPTER  V 


SUMMARY  AND  CONCLUSIONS 

Several  tests  of  the  finite  element  method's  ability 
to  handle  nonlinear  flow  problems  have  been  conducted.  The 
investigation  has  concentrated  on  spatial  discretization, 
consistent  versus  lumped  mass,  time  differencing,  and  com¬ 
putational  arrangement  of  the  equations.  When  combined,  the 
results  from  the  vorticity-stream  function,  penalty,  and 
shallow  water  equation  models  provide  consistent  answers 
about  the  application  of  the  finite  element  method. 

The  results  from  the  three  models  show  superior  per¬ 
formance  when  the  four-node  bilinear  element  is  used.  This 
conclusion  is  dependent  on  the  total  number  of  nodal  or  grid 
points  remaining  the  same.  For  example,  if  the  nine-node 
quadratic  elements  were  the  same  size  as  the  four-node  ele¬ 
ments  there  would  be  a  four-fold  increase  in  the  total  num¬ 
ber  of  grid  points  and  a  corresponding  increase  in  computer 
time;  however,  the  accuracy  would  be  greater  than  the  four- 
node  element  of  equal  size. 

Apparently,  the  key  to  element  accuracy  is  related  to 
the  influence  on  the  solution  at  a  given  node  by  the  surround- 
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ing  nodes.  Because  of  the  connectivity  of  the  elements,  the 
solution  at  a  node  is  affected  by  all  the  nodes  of  the  ele¬ 
ments  which  share  that  nodal  point.  For  example,  when  the  mesh 
consists  of  four-node  elements  each  nodal  point  belongs  to 
four  elements  so  that  the  solution  at  a  given  node  depends  on 
the  eight  surrounding  nodes.  This  is  true  at  all  interior 
nodes.  In  contrast,  a  mesh  made  up  of  nine-node  rectangular 
elements  works  differently.  The  corner  nodes  of  the  nine- 
node  element  are  part  of  four  elements  so  that  their  solution 
is  affected  by  the  twenty-four  surrounding  nodes.  The  mid¬ 
side  nodes  of  a  nine-node  element  are  part  of  two  elements  so 
that  fourteen  other  nodes  affect  the  solution.  And  finally, 
the  center  node  of  a  nine-node  element  is  only  influenced  by 
the  other  eight  nodes  in  the  same  element.  Other  elements 
have  similar  relationships.  In  Fig.  5,  it  is  evident  that  the 
triangular  element  can  result  in  different  nodal  relationships 
depending  on  the  element  orientation.  The  worst  case  occurs 
for  slant  5  where  the  influence  alternates  between  four  and 
eight  surrounding  nodes.  In  a  system  of  equations  which  al¬ 
lows  gravity  waves  this  inconsistent  influence  of  the  nodes  on 
each  other  could  enhance  the  development  of  gravity  waves. 

The  use  of  all  elements  of  one  size  or  shape  may  not 
always  be  possible,  especially  if  the  finite  element  method  is 
applied  to  boundary  layer  work  as  planned  by  Gresho,  et  al . 
(1978a)  or  to  models  where  it  is  desirable  to  have  a  nested 
fine  mesh  grid.  Further  research  is  required  in  order  to 
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determine  the  effects  of  variable  mesh  on  the  solutions.  How¬ 
ever,  this  study  indicated  that  the  rectangular  element 
should  be  the  primary  element  used  in  advective  flow  models 
with  the  triangular  element  employed  only  when  absolutely 
necessary  to  accomodate  an  irregular  boundary.  The  use  of 
higher  order  elements  is  neither  economically  nor  computa¬ 
tionally  desirable. 

This  study  substantiates  the  findings  of  Gresho,  Lee, 
and  Sani  (1977)  regarding  consistent  and  lumped  mass.  Study¬ 
ing  an  advective-diffusive  model  involving  dissipation  of 
pollutants  they  found  that  the  solution  suffers  when  mass 
lumping  is  applied.  Here  and  in  concert  with  the  related 
work  of  Sasaki  and  Chang  (1979) ,  it  is  clearly  shown  that 
regardless  of  the  time  differencing  scheme  used,  the  accuracy 
of  advective  calculations  is  severely  reduced  by  lumping  the 
mass.  Since  the  mass  is  inherently  lumped  in  finite  differ¬ 
ence  models  the  loss  of  accuracy  is  normally  improved  by  using 
a  finer  grid  resolution.  As  shown  by  Cullen  (1974a) ,  for  some 
problems  the  finite  element  method  with  consistent  mass  per¬ 
forms  as  well  as  a  finite  difference  model  with  four  times  the 
grid  resoltuion. 

Related  to  the  question  of  consistent  versus  lumped 
mass  is  the  problem  of  choosing  a  time  difference  scheme. 

Aside  from  the  problem  of  controlling  gravity  waves  it  is  ap¬ 
parent  that  when  consistent  mass  is  employed,  there  is  little 
difference  in  the  solution  as  long  as  the  time  integration  is 


91 


not  carried  for  a  long  time.  As  shown,  the  Crank-Nicholson 
scheme  allowed  long  term  time  integration  with  only  slight 
changes  in  the  total  energy  while  the  Leap-frog  scheme  al¬ 
lowed  faster  development  of  the  gravity  waves.  However,  with 
the  four-node  element.  Leap-frog  time  differencing  is  four 
times  faster  than  Crank-Nicholson.  Based  on  these  results, 
an  explicit  scheme  such  as  Leap-frog  appears  to  be  a  good 
compromise,  especially  if  it  is  combined  with  some  type  of 
time  filter.  Other  alternatives  include  the  use  of  the  semi- 
implicit  scheme  of  Kwizak  and  Robert  (1971)  as  planned  in  the 
model  under  development  in  Canada  (Staniforth  and  Mitchell, 
1978). 

The  final  topic  under  consideration  deals  with  the 
computational  arrangement  of  the  equations.  From  an  economic 
standpoint,  much  can  be  gained  by  rearrangement  of  the  equa¬ 
tions.  Computation  time  can  be  drastically  reduced  in  the 
finite  element  method  by  solving  the  equations  for  the  time 
rate  of  change  of  the  variables  rather  than  the  variables 
themselves.  By  using  additional  storage,  the  mass  matrix  can 
be  preprocessed  and  stored,  saving  additional  time.  If  the 
finite  element  method  is  applied  to  a  problem  with  uniform 
grid  spacing  then  it  is  possible  to  transform  it  into  a  quasi- 
finite  difference  model  thus  increasing  the  computer  effic¬ 
iency.  Such  a  transformation  has  been  made  by  Jesperson  (1974) , 
for  the  vorticity-stream  function  model  and  as  shown  here,  re¬ 
sults  in  a  four-fold  savings  of  computer  time.  Therefore,  even 


though  the  finite  elemerc  method  is  complicated  compared  to 
conventional  finite  difference  methods  it  can  be  made  very 
competitive  and  can  result  in  increased  accuracies. 


In  summary,  further  research  into  the  application  of 
the  finite  element  method  is  warranted.  In  particular,  the 
following  topics  should  be  investigated. 

1.  Because  of  the  potential  uses  in  meteorology,  a  de¬ 
tailed  examination  of  the  effects  of  grid  refinement 
should  be  made.  Several  options  such  as  reduction  of 
the  grid  through  the  prudent  use  of  triangular  ele¬ 
ments,  smooth  versus  abrupt  reductions  in  rectangular 
element  size,  and  moving  fine  mesh  grids  inside  coarse 
mesh  grids  should  be  studied. 

2.  The  use  of  variational  constraints  as  developed  by 
Sasaki  (1976)  should  be  investigated.  The  finite  el¬ 
ement  method  naturally  conserves  certain  properties 
depending  on  the  governing  equations;  however,  as 
demonstrated  by  Sasaki  using  finite  difference  methods, 
some  improvement  can  be  made  through  the  use  of  vari¬ 
ational  constraints. 

3.  Modelling  of  the  equations  on  the  globe  should  be 
studied.  Cullen  (1974b)  used  a  specialized  triangular 
mesh  for  this  work  but  did  not  perform  an  exhaustive 


study.  In  particular,  the  finite  element  method 
should  be  able  to  handle  the  changing  resolution 
caused  by  the  geometry  as  well  as  the  singularity 
points  at  the  north  and  south  poles. 

4 .  The  extension  of  the  model  to  the  vertical  dimen¬ 
sion  should  be  examined.  Staniforth  and  Daley  (1979) 
have  begun  development  of  a  baroclinic  model;  however, 
studies  should  be  made  as  to  the  required  resolution 
in  the  vertical  when  the  FEM  is  employed. 

5.  Investigations  should  be  made  into  the  modelling  of 
terrain  in  a  finite  element  model.  Preliminary  stud¬ 
ies  can  be  made  with  the  shallow  water  equations  in 
two  dimensions  by  placing  an  obstruction  in  the  chan¬ 
nel.  More  sophisticated  studies  should  then  be  made 
with  three  dimensional  flow  and  eventually  with  baro¬ 
clinic  models. 
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APPENDIX  A 


THE  FINITE  ELEMENT  METHOD 

Introduction 

The  purpose  of  this  appendix  is  to  provide  a  capsule 
description  of  some  of  the  pertinent  details  about  the  finite 
element  method.  Topics  discussed  are  the  Galerkin  method, 
discretization  of  the  domain,  and  development  of  interpola¬ 
tion  functions.  Additional  details  may  be  found  in  Segerlind 
(1976),  Strang  and  Fix  (1973),  Oden  and  Reddy  (1976),  and 
Zienkiewicz  (1977). 


Galerkin  Method 

The  Galerkin  method  is  a  variational  method  which  can 
be  applied  to  linear  as  well  as  nonlinear  problems  because 
there  is  no  requirement  that  the  functional  be  known. 

Suppose  the  following  governing  equation  is  specified 

-V2u  -  f  in  fl,  u  ■  0  on  30  (Al) 

The  solution,  u  (x,y)  is  approximated  by 

N 

u(x,y) ~  u  (x , y )  -  Z  u{N.(x,y)  (A2) 

6  i«l  1  1 
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where  are  approximating  functions  which  satisfy  some  con¬ 
tinuity  requirements  as  well  as  the  boundary  conditions.  If 
the  approximation  (A2)  is  substituted  into  the  governing 
equation  (Al)  the  result  is 


-V*u  -f  -  R(x,y) 

C 


(A3) 


where  R(X/y)  is  the  error  of  the  approximation.  The  Galerkin 
method  minimizes  this  error  as  follows 


JJ  R(x,y)  Ni(x,y)  dxdy  -  0, 

n 

or  JJ  (-V2ue~f  )Ni  (x,y)  dxdy  *  0. 


(A4 ) 
(A5) 


So,  the  error  is  made  orthogonal  to  the  set  of  approximating 
functions.  Substituting  u_  into  (AS)  gives 

l  JJ  (-V2  N.  (x,y))N  dxdy  u '■  JJ  f  N.  (x,y)dxdy.  (A6) 
j  a  J  1  3  a 

The  left  side  of  (A6)  can  be  transformed  using  the  product 
rule  of  differentiation  and  Green's  theorem  giving 


*  JJ  <_5ir  Tx  *  "jy  Ty 1  dxd*  uj  mJJ  1  »i<*-y)axdy 

+  f  ,ue"inx  *  ViY 
30 


(A7) 


.)  ds. 


where  (n  ,n  )  are  the  unit  normal  vectors  to  the  boundary 
x  y 

3(2.  Since  u  is  specified  (Al)  on  3(2  the  boundary  term  is 
zero  and  may  be  dropped.  Thus  (A7)  becomes 


[K]  {u>  -  {P}, 


(A8) 
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where 

tr  3N.  3N .  3N.  3N. 

Kij  =  //  ("?F  Tx  +  Tf  “iy }  dxdy' 

n 

=  jy  f  dxdy. 

ft* 

In  this  example,  the  Galerkin  method  was  applied  over  the 
domain  ft.  in  the  finite  element  method,  the  Galerkin  method 
is  applied  by  defining  the  approximation  (A2)  over  small 
elements  so  that  there  is  piecewise  continuity  within  the 
domain.  The  contributions  from  each  subdomain  are  summed  to 
approximate  the  global  domain. 


Discretization  of  the  Domain  and  Interpolation  Functions 


The  ways  of  discretizing  the  domain  are  only  limited 
by  the  imagination.  In  general,  the  governing  equations  or 
Galerkin  integrals  determine  the  minimum  order  of  the  approx¬ 
imation.  Higher  order  approximation  than  that  required  may 
be  used  and  in  some  cases  may  give  better  results.  The  geo¬ 
metry  of  the  domain  often  determines  the  shape  of  the  ele¬ 
ments  and  placement  of  the  nodes.  Fig.  A1  shows  a  domain  sub¬ 
divided  by  the  finite  difference  method.  Fig.  A2  shows  how 
the  same  domain  could  be  discretized  for  application  of  the 
finite  element  method.  Note  that  in  Fig.  Al,  there  are  sev¬ 
eral  grid  points  outside  the  domain  which  require  special 
handling.  For  the  finite  element  discretization.  Fig.  A2, 
two  elements  were  used.  The  domain  shown  could  represent 
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Fig.  Al.  Discretization  of  domain  for  finite  differences 


flow  in  a  channel  with  an  obstruction.  Fig.  A2  shows  some  of 
the  flexibility  of  the  finite  element  method  in  that  finer 
resolution  is  used  over  the  obstruction  and  coarser  resolu¬ 
tion  on  the  edges.  Obviously,  there  are  many  possibilities 
for  discretization  depending  on  the  problem  and  the  domain. 
Discretization  of  the  domain  leads  to  development  of  the  in¬ 
terpolation  functions. 


Fig.  A2. 


Discretization  of  domain  for  finite  element  method 


Finite  Element  Interpolation  Functions 


There  are  two  basic  requirements  for  interpolation 
functions;  completeness  and  compatibility.  Suppose  that  the 
element  equations  contain  derivatives  of  order  m.  The  com¬ 
patibility  requirement  means  that  at  the  interelement  bound¬ 
aries,  the  field  variable  and  its  partial  derivatives  up  to 
order  m-1  must  be  continuous.  The  completeness  requirement 
means  that  within  an  element,  all  the  constant  states  of  the 
field  variable  and  its  derivatives  up  to  order  m  must  be  re¬ 
presented  in  the  interpolation  function.  That  is,  the  poly¬ 
nomial  used  as  the  interpolation  function  must  contain  all 
the  terms  up  to  order  m. 

Suppose  we  have  a  system  of  equations  in  which  first 
order  derivatives  are  the  highest  order  appearing  in  the 
Galerkin  expressions.  The  compatibility  requirement  means 
that  the  minimum  order  for  the  interpolation  functions  is  one. 
Assume  that  the  global  domain  can  be  discretized  using  tri¬ 
angular  shapes.  We  can  choose  the  linear  approximation 

u(x,y)  *  ao  +  a^  +  a2y.  (A9) 

Note  that  if  a  node  is  placed  at  each  vertex  of  the  triangle 
there  are  three  unknowns  (aQ,  a^,  and  corresponding  to 
the  three  nodes.  This  polynomial  satisfies  the  compatibility 
and  completeness  requirements.  The  next  step  is  to  determine 
the  coefficients  aQ,  a^,  and  a2  for  the  arbitrary  element 
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rig.  A3.  Typical  triangular  element. 

shown  in  Fig.  A3. 

From  (A9 )  we  get 


or  {u}  =  [A]  {a>  , 

{*}  =  [A]'1  {u}, 


The  determinant  (det)  of  [A]  turns  out  to  be  two  times  the 
area  of  the  element.  At  any  point  (x,y)  in  the  element 


u  =  {l  x  y}  {a}  =  (1  x  v;  [A]  (u)  , 


(All) 


where 


u  =  ZN.u., 


N.  =  -i—  (r .  +  s.  X  +  t.  y) . 
i  det  11  i 


ri 55  Vk  -  Vj- 

si  =  -  v 

=  xk  -  ^  • 


i  t  j  *  k 


Other  Methods  of  Determining  Interpolation  Functions 


There  are  several  different  means  of  deriving  inter¬ 
polation  polynomials.  The  most  common  method  for  rectangular 
shaped  elements  is  to  multiply  the  one  dimensional  polynomial 
for  each  dimension.  That  is,  in  one  dimension  we  have 


u(x)  =  b  +  -b.x  or  u(v)  =  b_  +  by. 
o  1  *  i  o 


If  u  were  a  function  of  x  and  v  we  would  take  the  product 


u(x,y)  =  aQ  +  a^x  +  a2y  +  a^xy. 


Another  method  involves  the  use  of  Pascal's  criangie  [24] 


Fig.  A4 .  Pascal's  Triangle 


The  terms  connected  by  the  solid  line  represent  the  terms 
needed  for  a  quadratic  polynomial  to  be  used  on  a  triangular 
element.  There  are  six  terms  and  such  an  element  would  have 


i 

I 

six  nodes.  Note  that  the  dashed  line  encloses  nine  terms, 
these  nine  terms  would  be  used  for  a  nine-node  rectangular 
element.  The  location  of  the  terms  relative  to  the  dashed 
box  indicated  the  nodal  locations  in  the  rectangle. 

Local  Coordinates 

Previously  the  element  was  defined  in  terms  of 
the  global  coordinates  (x,y) ;  however,  it  is  often  conven¬ 
ient  to  define  the  element  in  some  local  coordinate  system 
and  later  transform  that  system  to  the  global  coordinates. 

This  is  especially  true  of  the  triangular  element  where  an 
area  coordinate  (Segerlind  1976)  is  convenient  (see  Fig.  A5) . 


Fig.  A5.  Area  coordinates  for  the  triangle 
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The  point  (£ ,  n)  in  Fig.  A5  divides  the  triangle  into 


three  areas  (L^,  ) .  We  can  define  the  shape  functions 

in  terms  of  these  areas.  Let  L  be  the  total  area  of  the 
element,  then 

Ni  =  Li/L  i=  1-3 

Lx  +  L2  +  L3  =  L 

Interpolation  Functions  for  the  Elements  Studied 


Linear  Triangle 


L1/L(2L1/L-1) 

L2/L(2L2/L-1) 

L3/L(2L3/L-1) 

4L1L2/L2 

4L,L-/L2 

4.  J 

4L1L3/L2 


APPENDIX  B 


LINEAR  STABILITY  ANALYSIS 

Introduction 

The  purpose  of  this  appendix  is  to  show  the  linear 
stability  analysis  for  the  time  differencing  schemes  used 
in  this  report.  Because  of  the  superior  results  achieved 
through  the  use  of  the  four-  node  linear  element,  the  stabil¬ 
ity  analyses  are  only  performed  for  that  element.  The  tech¬ 
nique  used  is  the  von  Neuman  method  as  discussed  by  Mesinger 
and  Arakawa  (1976)  . 


Stability  Analysis 

Consider  the  linear  advection  equation. 

If "  -  “o  H  -  vo  If  -  (B1> 

where  uq  and  vQ  are  positive  constants.  This  equation  des¬ 
cribes  the  advection  of  the  variable  $  at  a  constant  velocity 
c  given  by, 

c  *  (uQ2  +vol>*  (B2) 
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The  variable  <j>  is  assumed  to  be  of  the  form, 

<J>(x,y,t)  *  $(t)  exp  [i  (kx  +  ly)  ] ,  (B3) 

where  k  and  1  are  wave  numbers  in  the  x  and  y  directions, 
respectively. 

The  finite  element  method,  using  the  four -node  bi¬ 
linear  element  is  applied  to  (Bl)  as  discussed  in  Chapter  III 
and  Appendix  A.  By  expanding  the  matrices  in  the  resulting 
expression  a  finite  difference  form  of  the  finite  element 
equation  is  obtained.  The  resulting  expression  is 

yg-[16<Mj  ,m)  +  4  ($  ( j+l,m)+$  <  j-l,m)  +$  ( j  ,m+l)  +<Hj,m-l) 

+  $(j+l,m+l)  +£(j+l,m-l)  +<Mj-l,m+l)  +$ ( j-1 ,m-l) ] 

(B4 ) 

-u 

■  (<K  j-l*m+l)  -<Mj-l,m-l) )  +  (4>  ( j+l,m+l)  -<j>(  j+l,m-l) ) 

-v 

+4  ( <<>  ( j  , m+1 )  -0  ( j  , m-1 ) )  ]  [  (<f>  (j+l,m-l)  -<p  (j-l,m-l) ) 

+  (4i  ( j+l,m+l)-<j)  ( j-l,m-fl) )  +4  (<Mj+l,m)  -<j)  ( j-1  ,m)  ]  , 

where  Ax  and  Ay  are  the  grid  spacing  in  x  and  y,  respectively. 
This  expression  is  valid  at  any  interior  grid  point.  The 
terms  such  as  (j,m)  indicate  the  grid  point  where  x  =  jAx  and 
y  *  mAy. 

Leap-frog 

As  discussed  in  Appendix  A  and  Chapter  III,  the 
matrix  form  of  (B4)  is 
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[M]  (|}(n)  =  -  [K]  U)(n)  ,  (B5) 

where  [M]  is  the  mass  matrix  and  [K]  represents  the 

linear  advection  and  the  superscript  indicates  the  time  step. 
For  the  Leap-frog  scheme,  (B5)  becomes 

[M]  {q>}(n+1)  -  [M]  {4>}<n"1)  =  -2At  [K]  {$}  (n)  .  (B6) 


Eq.  (B3)  is  written  in  discrete  form: 
x  (n) 


(j,m)  =$  An  exp[i  (kjAx  +  ImAy) 1 , 


(B7) 


where  <J>Q  is  the  initial  amplitude  and  An  is  the  amplification 

A. 

factor  to  the  n  power.  Using  the  Eqs.  (B4,B6,  and  B7)  it 
can  be  shown  that 

An+1  exp[i(kjAx  +  ImAy)]  (2  +  cos  lAy)  (2  +  cos  kAx) 

-An_1  exp[i(kjAx  +  ImAy)]  i-  (2  +  cos  lAy)  (2  +  cos  kAx) 

jm  u  (B8) 

■  -  -g-  i  An  exp [i (kjAx  +  ImAy)]  [sin  lAy(2+cos  kAx)] 

v 

+  [sin  kAx (2  +  cos  lAy) ] }  . 

For  simplification,  assume  that  uo=vq  and  Ax  ■  Ay  ■  d.  Sub¬ 
stituting  (B2)  into  (B8)  and  simplifying  gives 

*’♦3/1  n  *  StS0*Ad>  -  1  -  »•  <B9> 

The  term  in  parenthesis  in  (B9)  has  a  maximum  value  of  2//3. 
Making  that  substitution  and  solving  for  X  gives 

(BIO) 


nr  cAt  ,  .  .  ,  cZAt2  .  ,  „  *s 
X  ■  -/6  -tj—  i  +  (-6  -jj +  1) 
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It  can  easily  be  seen  that 


| A | -  l  if  /E  £££  <  1.  (Bll) 

Therefore,  the  coupling  of  the  time  derivative  through  the 
consistent  mass  matrix  as  in  (B4)  results  in  a  more  restric¬ 
tive  time  step  than  is  found  in  finite  difference  applications 
or  with  lumped  mass.  A  similar  analysis  is  performed  by 
Cullen  (1974a) . 

If  the  left  side  of  (B4)  is  replaced  by  the  single 
term  <Mj,m)  then  the  stability  criteria  is 

/7  ^  <_  1.  (B12) 


"Theta-family"  of  Time  Approximations 


For  the  "theta-family"  of  time  approximations  dis 
cussed  in  Chapter  III,  (B5)  becomes 


lM]{<D(n+1)  -  [m]  (<D(n)  =  -eat  [k]  {$}  (n+1)-(i-e) At [K]{<i>}(n) . 

(B13) 

Using  Eqs.  (B4,  B7,  and  B13)  as  well  as  all  the  simplifica¬ 
tions  used  to  obtain  (BIO)  gives 


X 


_  1  ~  (1-9)  pi 

1  +  9pi 


(B14 ) 


where  p  =  /S’  . 

For  the  Crank-Nicholson  scheme,  0  *  1/2,  and 

I  A |  -  1. 
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Therefore,  the  Crank-Nicholson  scheme  is  always  stable.  For 
the  Backward  scheme,©  ■  1,  and 


M 


1 

1+p  * 


If  p  >  0,  this  scheme  is  always  stable  but  the  solution  is 
damped.  The  larger  p  is,  the  larger  the  damping  effect.  For 
the  Galerkin  scheme,  6=  2/3,  and 


1*1 


_  1  +  P/3 
"  r~4-  2p73 


This  scheme  is  also  stable  under  the  same  conditions  as  the 
Backward  scheme.  Damping  is  also  present  in  the  Galerkin 
scheme  but  the  damping  is  less  than  in  the  Backward  scheme. 
If  the  "theta-family"  of  approximations  is  analyzed  for 
lumped  mass,  the  only  difference  in  the  result  is  that  p  = 
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